Abstract
Understanding the microheterogeneity of tight sandstone is the basis of reservoir science, and quantitative characterization of the reservoir’s microheterogeneity is key to reservoir evaluation. In this study, an image-processing analysis method to study the heterogeneity of tight reservoirs is established. A modified Image J plugin is used to accurately identify the surface porosity of a thin casting sheet; the heterogeneity of the microscopic pores in a reservoir is then abstracted into discrete values of the surface porosity. A new parameter for quantitative characterization of the microscopic heterogeneity of tight sandstone, that is, the heterogeneity index
Introduction
Studies on the reservoir pore structure are the basis for formulating oil and gas field exploration and development plans, as well as reservoir evaluation (Noman and Archer 1987; Yang et al. 2013a; Yang et al., 2013b; Lim, 2017). With the increasing demand for energy in modern society, unconventional oil and gas resources are receiving increasing attention (Dou et al., 2021). With the exploration and development of unconventional oil and gas resources, the reservoir storage space has developed gradually to smaller levels and the study of the reservoir’s microscopic pore structure has become a top priority for current oil and gas resources development (Christopher et al., 2019; Zhang, 2018). At present, microscopic pore structure characterization methods for tight sandstone reservoirs are divided into qualitative and quantitative methods (Zhang et al., 2019). Qualitative characterization technology is primarily based on two-dimensional planar image observation scanning and three-dimensional spatial reservoir structure reconstruction simulations, which can reveal the preliminary shape, size, distribution, connectivity, oil–water occurrence state, and pore throat rock wettability (Li et al., 2018; Liu et al., 2004; Tian et al., 2014; Wang et al., 2017; Zhang, 2017; 2018).
To improve the characterization characteristics of the microscopic pore structure of reservoirs, many studies on microscopic pore quantitative characterization methods have been conducted. Characterization methods can be divided into three categories: (1) image analysis methods such as the observation and analysis of images obtained via optical microscopy and scanning electron microscopy; and (2) experimental measurement methods such as mercury intrusion, nuclear magnetic resonance, and gas adsorption. Accordingly, several scholars have conducted innovative research based on the abovementioned methods (Li et al., 2018; Sun, 2019; Yue and Liu, 1996).
Microscopic heterogeneity is an important factor that affects reservoir seepage and is fundamental to the microvalve pore structure (Shanley et al., 2004). More importantly, it affects oil recovery by affecting fluid seepage (Zhong et al., 2020). At present, the study on reservoir heterogeneity has made some progress in sequence stratigraphy, reservoir geological modeling methods, comprehensive geological methods, and various experimental methods. The geological cause analysis method is the method used most commonly for reservoir heterogeneity research; however, it lacks a quantitative analysis (Weger et al., 2009). In geostatistical methods, there are three commonly used parameters to quantitatively characterize the reservoir’s heterogeneity: the permeability coefficient of variation, the sudden increase coefficient, and the permeability (Tsakiroglou and Payatakes, 2000). However, these macroscopic heterogeneity characterization parameters not only ignore the change processes and the direction of the reservoir permeability but also have no upper and lower limits for their calculation, lack references, and are not convenient for quantitative evaluations. The biggest drawback of these parameters is that they are not backed by microscopic pore heterogeneity studies and ignore the effects of heterogeneities in the millimeter-, micron-, and nanoporosity on oil and gas enrichment.
In recent years, due to the development of fractal geometry, a new method to quantitatively describe the complexity of irregular heterogeneity in nature has emerged (Mahamud and Novo, 2008; Othman et al., 2010). This means that fractal theory can be used to classify and evaluate the heterogeneity of reservoir pore structures. Mandelbrot first proposed fractal theory and applied fractal geometry to compensate for the shortcomings of Euclidean geometry in the description of special structures. Subsequently, fractal theory has gradually been applied in the field of geosciences (Mandelbrot, 1982). Pfeifer and Avnir first discovered the fractal structure of reservoir pores via the molecular adsorption method and studied the adsorption fractal characteristics (Pfeifer et al., 1989; Pfeifer and Avnir, 1983; Pfeifer and Liu, 1997). Katz and Thompson (1988) studied the fractal properties of sandstone pore structures via scanning electron microscopy and high-pressure mercury intrusion experiments. Radlinski et al. (2004) studied the fractal characteristics and applications of sandstone pores. The fractal dimension has become an important parameter to quantitatively describe the surface roughness or irregular shape of microscopic pores (Krohn, 1988). Recently, fractal theory has been used to quantitatively characterize the aggregation behavior of wax crystals (Wang et al., 2020).
In view of the above, in this study, the face rate is obtained based on an image-processing method and a new parameter is proposed, namely, the heterogeneity index

Technology workflow.
Experiments
Experimental sample
The samples were collected from the third member of the Shahejie Formation in the Xiliu 10 fault block of the Huabei Oil Field. The thickness of the Es3 stratum is approximately 75 m. The Es3 stratum belongs to a delta plain meandering channel sand body in a water inlet system. The lithology comprises an interbed with unequal thicknesses of light gray fine sandstone and purple red mudstone. The lithology is primarily lithic feldspar fine sandstone; the quartz content in the debris is 46%–56%, the feldspar content is 27%–32%, and the rock fragments content is 14%–21%. The main components of the rock fragments are acidic and medium-based rock fragments, followed by sedimentary rock fragments. The grains are moderately weathered, with a sub angular and sub rounded shape, and sorting with a good, point-to-line contact between the particles. The cements were composed of calcite and shale. These cements include pore-cemented and contact-pore cemented. The pore space is not developed in the reservoir, including intergranular dissolved pores and intergranular pores with the diameter ranging from 0.01 to 0.05 mm. In this study, five tight sandstone samples from the Xiliu 10 fault block were selected (Figure 2).

Core samples from the Xiliu 10 fault block of the Huabei Oil Field.
Experimental methods
Low-temperature nitrogen adsorption experiments
Conventional physical properties of the five tight sandstone samples were measured to obtain their respective porosities and permeabilities. Low-temperature nitrogen adsorption experiments were performed, and experimental data for the specific surface, average pore diameter, and total pore volume were obtained.
The low-temperature liquid nitrogen adsorption experimental instrument (Figure 3) was a 3H-2000PS1/2 static capacity ratio surface and pore size-measuring instrument produced by Bethesda

3H-2000PS1 specific surface and aperture analyzer.
Instrument Technology (Beijing) Co., Ltd. First, the samples were weighed. The tight sandstone samples were ground into 40-mesh powders, weighed on a balance, and 1.5 g of the sample powder was placed in the sample tube. Next, each sample was purged and degassed. The temperature was set to 200°C, and the degassing time was 8 h. The third step was to test the samples. After the sample tube returned to room temperature, it was transferred to the test position, the software was opened to enter the test process interface, and the parameters were set and tested. The final step was to output the test results.
Acquisition of casting sheet images
Samples were collected in single polarization mode using a polarizing microscope. In the first step, 3 × 3 points were marked on the rock sample casting sheet. In the second step, at the selected positions, 3 × 3 face images were collected centered on the selected points. Five sets of samples were collected according to the abovementioned steps.
Therefore, a total of 81 × 5 face images were collected in this experiment. The 81 face images for each set of samples were independent of each other with no overlapping portions; these images were stitched into a complete sample casting sheet image (Figure 4).

Schematic of the sampling under the microscope.
Methodology for quantitative characterizing the microscopic heterogeneity
Surface porosity acquisition
A traditional binary image was obtained using an image-processing software that directly converted the RGB image into a grayscale image (Figure 5(a)) and then set a grayscale threshold to identify the pores (Sun et al., 2017; Zhang et al., 2009). However, in the process of the grayscale conversion calculation, due to the complexity of the rock and the color inhomogeneity of the cast rubber, and in the process of obtaining the binary image, the calculated values of the nonporous pixel points and the pore pixel points can be the same, resulting in recognition error (Lu et al., 2017; Figure 6(b)).

Traditional grayscale recognition and the RGB recognition matrix processing diagram (according to Lu et al., 2017) (a)Traditional grayscale recognition (b)Self-editing plug-in RGB recognition.

The pore extraction result: (a) the original optical microscope image; (b) the pore extraction map obtained using the traditional method; and (c) the pore extraction map obtained using the edited plugin.
To resolve this problem, we used a modified Image J plugin to accurately recognize the color images of the pores (Figure 5(b)). Image J is an open-source image-processing software based on Java. It was developed by the National Institutes of Health. It is widely used in medicine, biology, aviation, and other fields (Deng et al., 2012; Feng and Zhang, 2015; Yu et al., 2015; Zhang et al., 2017). As shown in Figure 5(a), the conventional image-processing process directly converts an RGB image into a grayscale image for image-processing software and then sets a threshold to identify the aperture. In contrast, this method uses threshold segmentation, and the calculated values of the nonporous pixel points and the pore pixel points are the same. Therefore, it is difficult to completely identify and separate the pores, causing recognition errors. In this study, the Image J plugin is modified such that the plugin uses the three RGB channels of the color-casting picture; then, according to the color of the casting glue, three thresholds are set and the pixels satisfying the three thresholds are recognized as the aperture pixel points and assigned a value of 0 (shown as black in Figure 6). Pixels that cannot be qualified arbitrarily are recognized as nonporous pixels and are assigned a value of 1 (shown as white in Figure 6). The three thresholds are limited at the same time, which greatly improves the recognition ability of the pore pixels, eliminates the case where the nonporous pixels have the same calculation value as the pore pixels, and accurately extracts the range of the blue pixels in the casting sheets.
In the core image used in this study, the cast body was blue. Via multiple tests, the RGB threshold of the blue aperture in the image was determined to be
The Image J software is open source with free plugins. We chose to use the voxel counter plugin. This plugin is used to capture threshold pixels within an image and to display the count and volume fraction (the ratio of the threshold pixels to all the pixels). For binary maps, black pixels are automatically assumed to be threshold points. It can be seen that the volume fraction obtained by running the plugin can be used to directly obtain the value of the surface porosity.
Figure 6(a) shows an image of one face collected under the microscope for sample Y-4, while Figure 6(b) shows the pore extraction map obtained using the conventional method. Figure 6(c) shows the result of the three threshold screenings of the edited software used in this study, that is, the pore extraction map. To enhance the contrast effect, the same area is denoted by a red box in the three panels in Figure 6, which is then enlarged and shown in Figure 7. After comparing the three panels in Figure 7, it can be seen that in the pore extraction map obtained using the conventional method, because the pores are recognized by only a single threshold, the calculated value of nonporous pixels are likely to appear the same as that of porous pixels. When the calculated values of the pixels are the same, aperture recognition error occurs and the image appears unclear with the contour of the aperture blurred, resulting in an inferior aperture range compared to the original image under the sample mirror (Figure 7(b)).

Local pore extraction results: (a) a partially enlarged view of the original image; (b) a partially enlarged image obtained using the conventional method; and (c) a partially enlarged image obtained using the edited plugin.
To solve this problem, we compiled a set of inserts suitable to accurately identify the pores of the core casting sheets. The plugin improves the method of converting RGB images directly into grayscale images in the traditional image-processing process by simultaneously setting different thresholds on the three RGB channels. This improvement optimizes the algorithm and eliminates pore recognition error in the binarization process. An image of a face processed using this method is characterized by a clean overall image, a clear outline of the pores, and a good correspondence with the contour range of the aperture of the original image (Figure 6(c)).
As can be seen in Figure 7(a), except for the pores of the blue cast body, the entire matrix is used. In the pore extraction map obtained using the traditional method, part of the matrix is also misidentified as pore, which makes the entire image unclear and blurs the pore contours (Figure 7(b)), causing the surface porosity of the subsequent calculated image to deviate significantly from the actual surface porosity. The aperture extraction map obtained using the modified Image J plugin essentially eliminates this recognition error. It can be seen in Figure 7(c) that the image is clean and the pore contour is clear; furthermore, the contour range of the pore has a good correspondence with the pore contour range of the original image. The surface porosity obtained using this method is very close to the actual surface porosity.
Microscopic heterogeneity index obtained using Image J software
At present, studies on reservoir heterogeneity have made some progress in sequence stratigraphy, reservoir geological modeling methods, comprehensive geological methods, and various experimental methods (Chen et al., 2017). The geological cause analysis method is used most commonly for reservoir heterogeneity research; however, it lacks a quantitative analysis. In geostatistical methods, three parameters are commonly used to quantitatively characterize reservoir heterogeneity: the permeability coefficient of variation, the sudden increase coefficient, and the permeability (Lu et al., 2017). However, these macroscopic heterogeneity characterization parameters ignore the change processes and the direction of the reservoir permeability and lack microscopic pore heterogeneity studies, ignoring the millimeter and micrometer scales and subsequently the influence of the heterogeneity grade and nanopores on oil and gas enrichment. To address the abovementioned deficiencies, this study seeks to quantitatively characterize the microscopic pore heterogeneity parameters. Based on an image-processing method and geostatistics, we propose a new parameter to quantitatively characterize the microscopic heterogeneity. According to geostatistics, each group of samples can be abstracted into a sample population and each surface porosity value extracted from a sample can be abstracted into a random variable, where several random variables form a sample. The statistical characteristics of the sample can be obtained, and the overall characteristics of the sample can then be deduced (Liang 2011). In this experiment, a stronger heterogeneity of the selected sample indicates a larger range of the surface porosity of the sample and a more discrete value. Therefore, the reservoir microscopic heterogeneity can be abstracted as the discreteness of the sample surface porosity values. Statistically, two parameters, the standard deviation and the variance, are often used to describe the degree of dispersion of a set of values. As the variance is the square of the data and is very different from the measured value, it is difficult to measure directly. Therefore, the standard deviation of the surface porosity value of the sample casting sheet was chosen as the microscopic heterogeneity index to quantitatively characterize the microscopic heterogeneity of the reservoir.
The sample size error is directly related to the accuracy of the inferred estimate. Overall, a larger sample size results in a smaller representative error of the statistical estimate. Generally, a sample number of greater than or equal to 30 is considered a large sample size. To better infer the estimated sample characteristics, the entire sample sheet was fully segmented under the same microscope conditions and 81 face images were collected to minimize the error in the statistical estimation. The standard deviation of the 81 face images in this experiment was calculated using equation (2), that is, the microscopic heterogeneity index is
Taking sample Y-4 as an example, the calculation process of the heterogeneous index Q is elaborated and its statistical significance is discussed.
First, 81 face images of the sample Y-4 were taken (Figure 8). Then, the sample face images were batch-converted into a pore extraction map using the modified Image J plugin (Figure 9), and the surface porosity was calculated (Table 1).

Sample Y-4 microscope below the hole image.

Pore extraction diagram of sample Y-4 following the Image J treatment.
Sample Y-4 surface porosity values.
Figure 10 shows a scatter plot drawn from the data in Table 2. It can be seen that the numerical distribution of the surface porosity of sample Y-4 is disordered and that there is no tendency to polymerize. The results can be fitted as follows: exponential fit,

Surface porosity scatter plot of the data for sample Y-4.
Heterogeneous index

Surface porosity numerical normal distribution of sample Y-4.
From the aforementioned determination, two conclusions can be drawn. First, the heterogeneity index
The microscopic heterogeneity index of the five sample sets was calculated using equation (2), and the specific values are shown in Table 2.
Results and discussion
Experimental results
Experimental results of nitrogen adsorption
Test results included the Brunauer–Emmett–Teller multipoint specific surface area, total pore volume of the Barrett–Joyner–Halenda desorption method, average pore diameter, and sample adsorption–desorption curve (Figure 12).

Nitrogen adsorption–desorption isotherms of the five tight sandstone samples.
The nitrogen adsorption–desorption curves of the tight sandstone samples measured in the experiment are shown in Figure 8. It can be seen in Figure 8 that the isothermal desorption curves of the samples show obvious desorption hysteresis and that separation occurs near the medium relative pressure (
According to the IUPAC classification schemes (Figure 13), the isotherms of samples Y-1, Y-2, Y-4, and Y-5 are classified as Type C, sample Y-4 is classified as Type B. This indicates that the pores of the Xiliu 10 fault block are primarily cylindrical with high openness and include crack-shaped pores.

Correlation between the hysteresis characteristics and the shape of the pores (according to IUPAC). (a),(b),(c), (d) and (e)Different types of Nitrogen adsorption-desorption isotherms curves.
Table 3 shows the routine experimental data of the samples in this experiment, including the reservoir’s physical properties and the specific surface, average pore diameter, and total pore volume of the samples determined in the low-temperature nitrogen adsorption experiment.
Test results of the conventional experiment for the tight sandstone Es3 in the Xiliu 10 fault block.
Fractal dimension based on nitrogen adsorption experiment
In the 1970s, the French scientist, Mandelbrot, developed fractal geometry to explain irregularities and instabilities; this type of geometry has a highly complex structure and has achieved remarkable results. It has subsequently been used widely in the field of geosciences (Hao et al., 2013; Zhang et al., 2016). Several researchers (Hao et al., 2013; Li et al., 2018; Sun et al., 2014; Zhang et al., 2013, 2016; Li et al. 2002; Ma and Lin 2012; Wen et al. 2007; Xie et al. 2014; Xiong et al. 2014; Zhao et al. 2014; Yang et al. 2015) have performed pore structure fractal dimension calculations for coal, shale, and low-permeability reservoirs and have proved that the fractal dimension can be applied to reservoirs for the quantitative characterization of the heterogeneity, pore throat distribution, and pore surface roughness (Angulo et al., 1992).
In this study, a low-temperature nitrogen adsorption experiment was used to determine the fractal dimension. The Frenkel–Halsey–Hill (FHH) model (equations (3) and (4)) created by Pfeifer was used to calculate the fractal dimension of the sample (Pfeifer and Colem, 1990). The specific algorithm is as follows. According to the experimental data of low-temperature liquid nitrogen adsorption, the adsorption–desorption isotherms of the five sample sets were plotted (Figure 12) and the hysteresis was determined. The pressure points at which the hysteresis ring separated and closed determine the pressure range of the capillary condensation zone. Within the pressure range of the capillary condensation zone, the x- and y-axes were used to fit the data via the least squares method (Figure 14) and the slope was obtained to yield the fractal dimension

(a) Fractal calculation results for the nitrogen desorption curves of Y-1; (b) Fractal calculation results for the nitrogen desorption curves of Y-2; (c) Fractal calculation results for the nitrogen desorption curves of Y-3; (d) Fractal calculation results for the nitrogen desorption curves of Y-4; (e) Fractal calculation results for the nitrogen desorption curves of Y-5.
Fractal dimension calculation results.
Here,
Fractal calculation results for the nitrogen desorption curves of the five sets of tight sandstone samples.
In Figure 14, the correlation coefficient between
Discussion
Test heterogeneity index Q
Many scholars have studied fractal theory and proved that this method can yield significant effects with respect to the characterization of the reservoir heterogeneity (Chen et al., 2017; Jiang et al., 2018; Li et al., 2006, 2018; Sun et al., 2017; Yang et al., 2014). In this study, the fractal dimensions of five sets of samples were calculated using a nitrogen adsorption experiment to test whether the heterogeneity index
Figure 15 shows the relationship between the microscopic reservoir heterogeneity index and the fractal dimension defined via image analysis. As can be seen in Figure 15, the heterogeneous index of the sample changes with the fractal dimension of the sample, as determined by the nitrogen adsorption experiment, and the two are positively correlated, with a correlation coefficient as high as 0.916. Therefore, the microscopic reservoir heterogeneity index defined via the image analysis method has a good characterization effect for the heterogeneity of the reservoir, pore throat distribution, and roughness of the pore surface. A larger

Relationship between the heterogeneous index Q and fractal dimension.
Geological significance of the heterogeneous index Q
Heterogeneous index Q and the reservoir’s physical properties
The reservoir’s physical properties are related to the nature of the rock and complexity of the pore structure, including the roughness, distribution, shape, and connectivity of the pore surface. Figure 16 shows the relationship between the heterogeneity index

Relationship between the heterogeneous index Q and reservoir’s physical properties.
Heterogeneous index Q and the pore structure parameters
To investigate the relationship between the sample heterogeneity index

Correlations between the heterogeneous index and the (a) specific surface area, (b) total pore volume, and (c) average pore diameter.
The microscopic heterogeneity index has a good correlation with the reservoir’s physical properties and pore parameters, indicating that the heterogeneity index can characterize and quantitatively evaluate the reservoir well. This classification provides a new direction for delineating the dominant horizon and guiding development and mining.
Conclusions
A modified Image J plugin was used to identify the pores of core casting sheets. The plugin discards the traditional method in which the RGB image is directly converted into a grayscale image and then threshold segmentation is performed in the image-processing process, resulting in the same nonporous pixel point and pore pixel point calculation value and leading to an aperture recognition error. Instead, three thresholds were set on the three channels of the RGB image and a segmentation calculation method was used to accurately identify the aperture. The voxel counter plugin was used to process the binary image, and the surface porosity of the core casting sheet directly indicated the result. Using this method to calculate the surface porosity of the casting sheet, a result could be obtained faster with higher quality. According to geostatistical methods, the sample was abstracted into a sample population and each time a face image was collected, it was abstracted into a random variable sample and all the collected surface porosity values were abstracted into a single sample. It was verified that the surface porosity values in the sample conformed to a normal distribution, indicating that the surface porosity values were statistically significant. The microscopic heterogeneity of the reservoir can be abstracted into the discreteness of the sample surface porosity values, and the standard deviation of the surface porosity values was defined as the heterogeneity index The fractal dimension can be used to quantitatively characterize the reservoir heterogeneity, pore throat distribution, and pore surface roughness. The heterogeneous index and fractal dimension of the sample were fitted and found to have a good positive correlation, indicating that the heterogeneous index has a certain geological significance. There is a good correlation between the heterogeneous index and the physical parameters, for example, the specific surface area, average pore diameter, and total pore volume of the sample, indicating that the heterogeneous index can characterize and quantitatively evaluate the reservoir well. This parameter provides a new basis for reservoir evaluation and reservoir classification and provides a new direction for delineating advantageous horizons and guiding development and mining.
Footnotes
Acknowledgements
We appreciate very much to Prof. Li Zhiping for his help.
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was funded by National Science and Technology Major Project (2017ZX05009005).
