Sandstone 3D compositional microstructure analysis with synchrotron-based multi-energy x-ray micro-CT

Three-dimensional (3D) characterization of materials microstructure is important in understanding their physical and chemical properties. In this paper, the compositional microstructure of a sandstone sample from Yaodian area of Yan’an in Erdors Basin is analyzed using synchrotron-based multi-energy x-ray micro-CT data. The x-ray CT data of the sample are acquired at beam energies 25 KeV, 35 KeV and 45 KeV. The strip artifact in the sinogram images are corrected with a least square fitting method. The corrected sinogram image data are reconstructed with the X-tract software, and the reconstructed images under different energy spectra are calibrated. Finally, the multi-energy least square segmentation method is used to characterize the 3D microstructure of sandstone sample with the Data-Constrained Modeling (DCM) software. The results show that the sample has 71.8% of quartz and albite which is displayed as blue, calcite (22.3%) is displayed as red, pyrite (0.6%) is displayed as yellow and porosity (5.3%) is displayed as white.


Introduction
Quantitative knowledge of material microstructures is important in various science and engineering applications. For example, the microstructure of metallic materials can determine the macroscopic properties of metals; the microstructure of sandstone reservoir, especially the pore structure, plays a decisive role in oil and gas reservoir productivity [1,2].
Recently, spectral CT based on a photon counting detector (PCD) has attracted an increasing attention. This imaging mode can obtain the projection data of multiple energy channels simultaneously by setting a threshold value. Therefore, it can obtain the internal three-dimensional (3D) structure of the sample and improve the identification accuracy of the material components. Now it has been applied in the field of material analysis and biomedical research [3,4]. However, spectral CT based on PCD is still in the initial stage of research and development, imaging speed is slow. It is not commercially ready for industrial applications. While the synchrotron-based x-ray micro-CT, combining the advantages of synchrotron radiation x-ray and nondestructive examination, can obtain the monochromatic x-ray projection data of the sample, eliminate the hardening artifacts, avoid the CT value drift of the conventional CT and get the accurate CT value. The technique has been widely used in many fields such as biomedicine and materials science [5,6].
Multi-energy x-ray micro-CT based on synchrotron radiation can obtain multiple single-energy synchrotron radiation CT data which can make full use of energy spectrum information, eliminate beamhardening artifacts and can effectively improve the identification accuracy of components. Therefore, in this paper, 3D microstructure characterization of a sandstone sample was studied based on synchrotron radiation multi-energy x-ray micro CT.
In the study of microstructure distribution of a sandstone sample, the quality of the reconstructed CT images is essential. Many factors could affect the quality of reconstructed images, such as the stability of synchrotron Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. radiation light source, the linearity and pixel homogeneity of the detector, the mechanical accuracy of the rotation stage, and the accuracy of reconstructed algorithms. With the advancement of instrumentation and sophisticated reconstruction algorithm, some of these problems could be practically compensated. The effect of the x-ray light source intensity variation is becoming more evident. This is partly due to the long distance between the source and detector and the relative narrowness of the light beam. A small angular shift could translate to significant intensity variation at the detector.
The errors caused by light source such as intensity unevenness, super-saturation and penetration have been studied for synchrotron radiation CT technique, and its three basic error forms are given [6]. In this paper, the problem of the sinogram artifacts caused by inconsistency of light source intensity in space and time is analyzed and the correction method is proposed.
After the CT images have been reconstructed using the corrected sinograms, the characterization of the compositional distributions will be conducted. The image segmentation techniques have been widely used for determining the compositional microstructure for x-ray CT. The existing image segmentation techniques can be classified into two major approaches: a statistical approach and a structural approach. In the statistical approaches, the least square method is a simple and practical technique. So, in this paper, the least square segmentation method was used to obtain the 3D compositional microstructure of a sandstone sample. For synchrotron-based multi-energy x-ray micro-CT, the conventional image segmentation technique would give multiple different compositional maps, one for each energy. An elevated segmentation accuracy would be anticipated if the multiple reconstructed CT data sets would be combined. Therefore, a multi-energy least square segmentation technique was utilized to obtain a 3D compositional map by combining multiple reconstructed CT data sets in this paper.
The rest of this paper is organized as follows. In section 2, a sandstone sample and the acquisition of its multienergy x-ray CT data are introduced. In section 3, firstly, we analyze the sinogram artifacts caused by the inconsistency of synchrotron radiation source and present a correction method of this artifacts. Then the corrected CT sinogram images are reconstructed, and the reconstructed images under different energy spectra are calibrated. A multi-energy least squared segmentation technique is utilized to obtain the 3D microstructure of this sandstone sample by combining the two reconstructed CT data sets in section 4. The conclude of the paper is presented in the last section.

Sandstone sample and acquisition of measurement data
A sandstone sample with a diameter and a height of 3 mm and 20 mm was extracted from a core which has a diameter and a height of 100 mm and 200 mm, as shown in figure 1. The core was drilled from Yaodian area of Yan'an in Erdors Basin, where the geographic structure and layer belong to the eastern area of Shanbei Slope of Member 6 in Yanchang Formation of Triassic in Erdors Basin. Erdors basin is an important energy base in China with the features of wide distribution, multiple oil-containing-layers, high reservoir thickness as well as low porosity and permeability. In this area, the main oil reservoir is the Triassic Yanchang Formation, in which intergranular pore and dissolved pore are the two typical types of reservoir Pore that occupy 42.47% and 56.37% in the total of pore respectively. The main minerals of the sample are quartz, albite, calcite and pyrite.
Multi-energy x-ray micro-CT experiments were carried out on the BL13W beamline at the Shanghai Synchrotron Radiation Facility (SSRF). SSRF of China is a world advanced third generation synchrotron light source and the energy of the electron beam in the SSRF electron storage ring is 3.5 GeV electron volts. It ranks fourth in the world, behind the world's only three high-energy light source, the United States, Japan, Europe. SSRF was officially opened to research users in 2010.
The selected beam energies are 25 KeV, 35 KeV and 45 KeV. The exposure time for each image was 1 s, 1.5 s and 8 s at beam energies 25 KeV, 35 KeV and 45 KeV, respectively. The distance between the sample and the CCD detector was set to 40 mm. The combination of the CCD detector (with a native pixel size of 7.4 μm) with an optical lens (2×) provides an effective pixel size of 3.7 μm. For each CT scan, 720 transmission images (terminate x-ray intensity) were recorded during an 180°rotation with a fixed angular displacement of 0.25°. In the experiments, dark-field images and flat-field images (incident x-ray intensity) were also recorded. Every two flat-field images were recorded after every 40 transmission images. At the end of the experiments, dark-view fields were recorded. The size of each image is 1377×721 pixels. Figure 2 gives flat-field images 1, 17, 35 at energy 35 keV and figure 3 gives transmission images 40, 360, 720 at energy 35 keV.

Artifacts analysis and correction
Beer-Lambert Law tells us )represents the linear absorption coefficient of an object to be measured, I 0 is incident light intensity, I is the terminated light intensity and p is projection.
The ideal intensity of the incident light during the experiment is very difficult to measure. For purpose of CT reconstruction, I 0 is replaced with the flat-field images. The sinogram image can be obtained from the flat-field images and transmission images as: Figure 4 gives 400 th sinogram images at energy 25 keV, 35 keV and 45 keV. As can be seen from figure 4, there is almost no artifact when the energy is 25 keV. While the strip artifacts appear in the sinograms when the energy is 35 keV and 45 keV, and the higher the energy, the more serious the strip artifacts. The reason is as following.
In formula (2), incident light intensity I 0 is replaced by flat-field image, in order to make flat-field image and I 0 as close as possible, multiple flat-field images were recorded. In the experiments, two flat-field images were   recorded after every 40 transmission images or every 10 degrees of sample rotation. When the beam energy is low, the exposure time is short and the illumination variation is low. That is, the difference between flat-field image and the ideal data I 0 is small. So we can use flat-field image to replace I 0 without significant decrease in accuracy. But when the energy becomes larger, the exposure time become longer and the difference between the experimental and ideal data is bigger. In this case, the replacement may result in artifacts appearing in the sinogram images.
The difference between two flat-filed images reflects that the intensity of light source at different time is different. For one flat-field image, the value in the horizontal direction reflects the intensity of the light source received in a certain layer, and the value in the vertical direction reflects the intensity of the light source received in different layers. Figure 2 gives the flat-field images collected at different times. From figure 2, it can be seen that the intensity of the light source is different in different acquisition time and different space position. And its intensity changes are little in the horizontal direction, while changes are greatly in the vertical direction. So, we assume that the intensity of the light source is uniform horizontally in this paper. The focus of the correction is for intensity variation of light source in different layers and at different time. Figure 5 shows the row average density curves of the 1, 5, 9, 13, 17, 21, 25, 29, 33-th flat-field images at the three different energies, respectively. The horizontal axis represents the number of layers, that is, the longitudinal position of the sample, and the vertical axis represents the average intensity of the light source in a layer. The different curves represent the light source intensity of the flat-field images collected at different time. It is obvious from figure 5 that the light intensity received by different layers is different, and when the energy become large, the interval of adjacent curves increases, that is, the intensity difference of adjacent flat-field images becomes large. Similar to the flat-field image, the spatial and temporal instability of the light source intensity appears in the transmission images, as shown in figure 3. In general, when the energy is high, the exposure time is long, resulting in a greater variation in the intensity of the light source, which ultimately results in more artifacts in the sinogram images, as shown in figure 4.
The spatial and temporal inconsistency of the synchrotron light source will bring errors to the sinogram images, which will affect the accuracy of the material absorption coefficient in the reconstructed images and then affect the accuracy of the sample microstructure prediction.
A natural way to solve this problem is to standardize the intensity of the light source. Firstly, the intensity information of the light source is obtained from the flat-field images and the transmission images. Here, a digital image can be expressed as matrix A=(a ij ) I×J . Assuming that the intensity of light source is uniform in the [¯¯¯]  represents the light source intensity in the vertical direction. For transmission images, the row average density curve not only reflects the intensity information of the light source, but also reflects the information of the sample. In order to get the trend of the light source in the transmission images, the row average density curve near both the left-and right-hand sides of the transmission images which mark as A sidē was plotted. The image areas were selected such that they do not overlap with the sample projection image at any angle of rotation. In order to remove the noise and other additional factors on the light source intensity, the least square method is used to fit the curve of the average intensity of the light source for each flat-field image and transmission image. Then, the fitting curve is used as a correction standard to standardize each original image. The pseudo-code of the correction algorithm based on least square fitting is shown in table 1. Figure 6(a) gives the row average density curves for the flat-field images 1, 17, and 35. And figure 6(b) shows the fitting results for the three curves in figure 6(a). The corrected flat-field images and transmission images are shown as  After all flat-field and transmission images were corrected by the same methods, the sinogram images were again obtained using formula (2), as shown in figure 9. As can be seen from figure 9, the strip artifacts in the original sinograms at 35 KeV and 45 KeV are well removed. Because there is no artifact in the sinogram images when the energy is 25 KeV, only the flat-field and transmission images at 35 KeV and 45 KeV are standardized. Table 1. Pseudo-code for standardized algorithm of light source intensity.
for k=1:K (K is the total number of flat-field or transmission images) Computing the row average density Ā / A sidē Fitting Ā / A sidē using curve fitting method of least squares polynomial, obtaining the fitting curve Ã / A sidẽ Using Ã / A sidẽ standardize each original image k=k+1 End

CT image reconstruction and calibration
One advantage of the synchrotron light source is that it can be approximated as parallel beam. After appropriate preprocessing, a set of multi-slices parallel-beam projection are obtained. Among all the available CT reconstruction algorithm, the most computationally efficient and widely used one is the FBP-type algorithm. So a standard parallel-beam FBP algorithm is used to reconstruct each sinogram image. The parallel beam FBP reconstruction of all sinograms in this paper are accomplished using the X-tract software.
In order to ensure the accuracy of microstructure reconstruction, the edge part of the reconstructed image was removed, and the middle part of 600×600×700 part was remained. Then the reconstructed CT images were fully aligned with one another in three dimensions. The alignment is achieved by translation and rotation of x-ray CT images to maximize the cross-correlations between images at different beam energies. Figure 10 shows the 360-th layer final CT reconstruction images after alignment at 25 KeV and 35 KeV, the size is 600×600 pixels.

Multi-energy least square segmentation method 4.1. Method
The goal of this paper is to gain the compositions microstructure distribution of a sandstone sample. While the multiple spectral CT has gained considerable interests because of its energy-resolution in identifying and discrimination materials. For multi-energy X-ray CT, if the traditional least square segmentation method is used, the multiple microstructure distributions of components under multiple energies can be obtained for the same sample. The precision of segmentation will be improved if the microstructure distribution of components under multiple energies are combined. So in this paper, a multi-energy least square segmentation technique was utilized to obtain a compositional map. This approach is similar to the dual energy x-ray CT techniques, which utilise the fact that the different components have different energy-dependence behaviour in their x-ray absorption properties [7].
The model is defined based on a voxel. Let N is the total number of voxels, M is the total number of compositions in the system, a m denotes the m-th composition, and n a   n a n a n a a 0 corresponds to pore. When we get the value n a n m ( )(m=0, 1,K, M, n=1, 2,K, N) for each voxel, we get the microstructure distribution of the sample.
In the multi-energy least square segmentation technique proposed in this paper, it is assumed that there is only one material in each voxel, i.e. n a = 0 m a n a m a n a m a n a m m a n a m a n a m a n a m m a n a m a n a m a n a m

Experiment results and analysis
In order to obtain the main compositions in the sandstone sample, scanning electron microscope (SEM) was used. The results showed that quartz, albite (Na-Feldspar), calcite and pyrite were the main compositions, as shown in figure 11. It was found to be difficult to discrete quartz and albite, because the x-ray absorption properties of them are very similar. So they are taken as a kind of mixture with ratio of 1:1 in this paper and the mixed density of the mixture is 2.635 g cm −3 . The density of pyrite and calcite are 4.9 g cm −3 and 2.71 g cm −3 , respectively. The chemical formula, density and attenuation coefficient of the three components are shown in table 2.
After alignment, 2 CT reconstruction images with size 600×600×700 at 25 keV and 35 keV were used for multi-energy least square segmentation. The objective function (5) is solved under the conditions (4) and (6),  which is calculated by the Data-Constrained-Modelling (DCM) software [8,9]. The software can be downloaded from the website http://research.csiro.au/dcm. Figure 12 shows (a) the composition distribution at 360-th layer and (b) 3D microstructure distribution of the sandstone sample. In figure 12, blue denotes for mixed composition by quartz and albite, red denotes for calcite, yellow denotes for pyrite, and white denotes for void.
The overall volume fractions of the compositions in the sample were shown in table 3. As can be seen from figure 12 and table 3, the mixture of quartz and albite accounts for 71.8% by volume, calcite accounts for 22.3%, pyrite accounts for 0.6% and porosity accounts for 5.3%. The predicted porosity value of 5.3% is in fair agreement with measured value of (6±0.1)%. The lower porosity value may result from the assumption that a unique material is associated with a voxel. Some pores which are smaller than the CT resolution may not have been accurately accounted for.

Conclusion and discussion
In this paper, the microstructural characterization of a sandstone sample mainly consisting of quartz, albite, calcite and pyrite was studied by multi-energy micro-CT technique based on synchrotron radiation. In order to solve the problem that the quality of CT images is influenced by the instability of the synchrotron radiation light intensity, a normalization method was proposed to compensate such variations. This method firstly fits the intensity graph of the synchrotron radiation light using a least-square polynomial curve fitting method. And then the fitted curve was used to normalize the illumination variation on the flat-field and transmission images. These corrected sinogram images were used to carry out CT image reconstruction. In order to make the reconstructed images under different energy have maximum correlation, the reconstructed images under different energy are aligned. Finally, a multi-energy least-square segmentation approach has been used to model the microscopic compositional distributions. The predicted porosity value is in fair agreement with measured value. In this paper, it has been assumed that each voxel contains only one mineral phase. This assumption should be reasonably for sandstones. However, it may not always be adequate for other rock types. Work is in progress for other reservoir rock types without such assumption.
Microscopic reservoir pore texture refers to the shape, size and distribution of pores and throats and their interconnectivity in reservoir rocks, which strongly influence the petrophysical properties of reservoir rocks, especially in those tight sandstone reservoirs with nanoscale pore-throats (pores and throats), microscopic porethroat texture controls their porosity and permeability mainly. Therefore, enhancing the accuracy of picture scanned by x-ray CT has been the key to the research of tight reservoir. The synchrotron-based multi-energy x-ray micro-CT proposed in this paper can improve the accuracy of data scanned and provide a more accurate