Superpixel Segmentation Based on Spatially Constrained Subspace Clustering

Superpixel segmentation aims at dividing the input image into some representative regions containing pixels with similar and consistent intrinsic properties, without any prior knowledge about the shape and size of each superpixel. In this paper, to alleviate the limitation of superpixel segmentation applied in practical industrial tasks that detailed boundaries are difficult to be kept, we regard each representative region with independent semantic information as a subspace, and correspondingly formulate superpixel segmentation as a subspace clustering problem to preserve more detailed content boundaries. We show that a simple integration of superpixel segmentation with the conventional subspace clustering does not effectively work due to the spatial correlation of the pixels within a superpixel, which may lead to boundary confusion and segmentation error when the correlation is ignored. Consequently, we devise a spatial regularization and propose a novel convex locality-constrained subspace clustering model that is able to constrain the spatial adjacent pixels with similar attributes to be clustered into a superpixel and generate the content-aware superpixels with more detailed boundaries. Finally, the proposed model is solved by an efficient alternating direction method of multipliers (ADMM) solver. Experiments on different standard datasets demonstrate that the proposed method achieves superior performance both quantitatively and qualitatively compared with some state-of-the-art methods.

I. INTRODUCTION
Superpixels are alternative primitives to represent the original pixels in an image based on certain intrinsic properties, such as coherent spatial relationship, and similar appearance. Instead of utilizing pixel entities directly, superpixels can provide semantic representation and reduce the dimension of primitives [1], [2], which have been widely applied in practical and industrial tasks of computer vision, such as saliency detection [3]- [20], feature selection [21], object recognition [22], background subtraction [23], image enhancement and reconstruction [24]- [26]. In recent years, numerous methods with diverse problem formulations have been proposed for superpixel segmentation [27]- [31] and achieved excellent performance in theoretical study. However, one obvious weakness of the existing superpixel methods is the limitation to adapt to the local image details, such as the content boundaries and object contours. And in many cases of practical problems and industrial tasks [3]- [5], [22], the local details are the key information, which should not be ignored. For most methods, the detailed boundaries of image content are still hard to be well preserved, due to the lack of the prior knowledge about the shape and size of the superpixels in an image. These methods have to increase the number of superpixels and reduce the superpixel size to achieve the detailed boundaries. This may usually lead to large data redundancy in the sparse area, which is unwillingness for practical tasks. How to balance the detailed information and superpixel number is a challenging problem for superpixel segmentation to be applied in industrial tasks.
To alleviate this problem, we suppose each superpixel with independent information comes from one distinct subspace, and formulate the superpixel segmentation as a subspace clustering problem [32]- [34]. Subspace clustering is an effective way to segment data drawn from different dimensions or subspaces, and has achieved impressive performance in some real-world tasks, such as face clustering [34] and foreground extraction [35]. Moreover, it also has been adopted in industrial applications, e.g., hyperspectral image classification based on superpixel-level subspace learning [36]- [38]. For superpixel segmentation, each superpixel is regarded as a subspace containing pixels with similar intrinsic attributes, such as similar texture, boundary, color and coherent spatial information, while pixels in different subspace have different properties. Then the solution to the subspace clustering corresponds to the segmentation of superpixels. By this way, the pixels in the unit containing independent semantic information are more probably segmented into one superpixel, no matter how small the unit area is. At the same time, the size of the superpixels in sparse area are relatively large according to the similar attributes of pixels. Consequently, the detailed content boundaries and the superpixel number can be better balanced, i.e., using less superpixels to preserve more detailed information.
Nevertheless, it is still a challenging work by inducing subspace clustering theory to superpixel segmentation directly. Different from the general sample data in clustering, pixels in an image have certain spatial correlation that only the similar pixels in adjacency tend to be clustered as one superpixel. Directly applying the conventional subspace clustering models to superpixel segmentation task can not excellently work since the spatial correlation of pixels is not taken into consideration. This may cause the confusion and segmentation error of content boundaries. In other words, we should consider the spatial constraints when using the subspace clustering for superpixel segmentation task. To this end, we devise a spatial regularization term to emphasize the spatial correlation, and propose a spatially constrained subspace clustering based superpixel segmentation model to generate superpixels with more accurate and detailed boundaries, which is more appropriate for practical and industrial tasks.
In this paper, considering the limitation of superpixel segmentation in the application of practical industrial tasks, we propose a novel spatially constrained subspace clustering based superpixel segmentation model to produce more adaptive superpixels. In summary, the contributions of this paper include: • Considering the challenge of pursuing detailed boundaries, we formulate the superpixel segmentation task as a subspace clustering problem to generate content-aware supeprixels with more detailed boundaries and fewer primitives. To the best of our knowledge, this is the first attempt to induce subspace clustering theory to solve the superpixel segmentation problem, and effectively alleviate the challenge in superpixel segmentation. • To enforce the spatial adjacent pixels with similar attributes to be clustered into a superpixel corresponding to a subspace, we design a spatial regularization term in the loss function to constrain the local adjacency when clustering pixels. The ablation study has verified the effectiveness of the spatial regularization term. • The extensive experiments and discussions on the widely used benchmarks and datasets demonstrate that the proposed method achieves superior performance both quantitatively and qualitatively. The paper is organized as follows. The related works are presented in Section II. In Section III, we propose a novel superpixel segmentation model. The comparative experimental results and analyses are presented in Section IV. Finally, Section V concludes this paper.
II. RELATED WORK Various superpixel segmentation algorithms have been proposed for different motivations and tasks. They can be classified into two categories: unsupervised approach and supervised approach.

A. Unsupervised Superpixel Segmentation Methods
Unsupervised approach mainly includes two categories: graph-based methods and clustering-based methods.
Normalized cut (Ncut) algorithm [39] is one of the representative graph-based methods, which first groups pixels into several large areas by normalized cut, and then adopts Kmeans to segment the large areas into small units. However, it suffers from high computational complexity cost and low boundary adherence, which limit its utilization. Considering the image topological structure, topology preserved regular superpixel (TPS) [2] algorithm was put out to produce relatively regular superpixel. However, it suffers from low adherence to boundaries due to the demand of regular superpixel shape. In contrast to TPS, lazy random walk (LRW) [27] algorithm achieves high object boundary adherence by iteratively optimizing an energy function based on texture measurement. Recently, approximately structural superpixels (ASS) [31] was proposed to generate approximately structural superpixels by an asymmetrically square-wise superpixel segmentation way, which can largely reduce data amount as well as preserve image content boundaries.
Turbopixel [40] is one of the early clustering-based methods, which adopts local image gradient and geometric flow to generate regular superpixel, but it suffers from the time consuming problem in practice. Simple linear iterative clustering (SLIC) [28] is widely used due to its high efficiency. SLIC adopts K-means clustering approach and changes the search area into a local area rather than the global searching in Kmeans, which largely reduces the computation complexity. As an extension of SLIC, Fast linear iterative clustering (FLIC) [30] is an active search method which emphasizes the neighboring continuity. Linear spectral clustering (LSC) [29] generates superpixels by mapping pixels into a high-dimensional space through a kernel function. Recently, Bayesian Adaptive Superpixel Segmentation (BASS) [41] is proposed to generate adaptive superpixels by a Bayesian mixture model. However, the number of superpixel cannot be predefined, which is not convenient for certain tasks that demand to initialize the number manually.

B. Supervised Superpixel Segmentation Methods
Supervised superpixel segmentation algorithms are mainly based on deep-learning based approaches. Recently, deep learning based methods for superpixel segmentation have been proposed to incorporate deep learning technique with traditional segmentation. There are three representative deep learning based superpixel segmentation methods, Segmentationaware Loss(SEAL) [42], Superpixel Sampling Networks(SSN) [43], and Superpixel Segmentation with Fully Convolutional Networks (FCN) [44].
SEAL adopts deep neural network and designs a segmentation-aware affinity learning method to keep weak object boundary in superpixel segmentation. SSN presents a differentiable superpixel sampling model which can be integrated into end-to-end trainable networks. FCN algorithm adopts a simple fully connected network to predict superpixel on regular image grids. Then, on the basis of the predicted

Global regularization Spatial regularization
Error Matrix ... × = + Fig. 1. The framework of our proposed method. Given an input image, we first construct the input data as a locality-constrained subspace clustering model with spatial regularization, and then work out the coding matrix by the proposed approach. Finally, the superpixels are generated via clustering the affinity graph of the coding matrix.  [28], FLIC [30], and LSC [29], respectively. (e) shows the result based on a naive subspace clustering method. We also show the result using our spatially constrained subspace clustering method in (f).
superpixel, it further develops a down-sampling/up-sampling scheme for deep networks to generate high-resolution output for intensive prediction tasks. By this way, the stereo matching network architecture can predict both superpixel and parallax.

III. PROPOSED METHOD
In this section, we first explore a naive subspace clustering model in superpixel segmentation, and analyze the irrationality of the simple integration. Then, we propose a novel localityconstrained subspace clustering model with spatial regularization for image superpixel segmentation. An overview of the proposed method is illustrated in Fig. 1.

A. Superpixels via Naive Subspace Clustering
Superpixels aims at representing the original image by some units containing pixels with similar properties, such as close appearance and spatial correlation. Subspace clustering finds diverse subspaces with different dimensions and bases, and clusters each data point in one subspace [34]. Each superpixel can be regarded as a subspace containing pixels with similar intrinsic attributes, such as similar texture, boundary, color and coherent spatial information, while pixels in different subspace have different properties. Then, superpixel segmentation problem can adopt subspace clustering theory to find the segmentation of these subspaces.
A naive way to solve this problem is to directly apply off-the-shelf subspace clustering algorithms (e.g., [34]) to superpixel segmentation by constructing an affinity matrix by learning effective representation coding. For efficiency, we use the extracted features from the raw pixel units rather than the raw pixel features as input. Here, a self-expressive model [34] is adopted to learn the affinity matrix for efficiency, and the data is explained by itself: a data point can be represented as a linear combination of other data points in the same subspace, i.e., X = XZ + E, where X = [x 1 , x 2 , · · · , x n ] is the feature matrix, Z is the coding matrix, and E is the error matrix. Consequently, we first formulate the superpixel segmentation problem as: where λ is a trade-off parameter, and 0 n ∈ R n×1 is an all zeros vector. The constraint diag(Z) = 0 n is used to eliminate the trivial solution, i.e., Z equals to an identity matrix. The term Z 0 is a regularization to enforce the sparsity, representing the global subspace structure in X. 0norm means to find the sparsest solution. However, minimizing (1) is extremely difficult since it corresponds to an NP-hard problem. To find a nontrivial sparse representation efficiently, the 0 -norm is usually relaxed to its convex surrogate, i.e., the 1 -norm [34].
Nevertheless, directly applying the conventional subspace clustering models to superpixel segmentation task does not work well. Fig. 2 shows the segmentation results of some stateof-the-art superpixel methods. From Fig. 2(b)-(d), we can see that the spots in the butterfly wing are almost unrecognizable. In Fig. 2(e), we adopt the subspace clustering model in (1) to generate superpixels, and the boundaries of the spots are also not achieved. The reason is that different from the general sample data in clustering, image pixels have inherent spatial correlations, which is not taken in the consideration of the conventional models. Therefore, we design a locality constraint to pursue the spatial correlation.

B. Superpixels via Locality-Constrained Subspace Clustering
Considering the spatial correlation of pixels in an image that only the similar pixels in adjacency tend to be clustered as one superpixel, we devise a spatial regularization term to 3. A simple case of adjacent matrix. If two units are neighbors, the corresponding element is labeled as "1" in the adjacent matrix; otherwise, the element is labeled as "0". enforce the adjacent correlation. Then, we propose a localityconstrained subspace clustering model to generate contentaware superpixels with more detailed boundaries.
Spatial Regularization: Considering a data set X = [x 1 , x 2 , · · · , x n ] drawn from input image units, for image superpixel segmentation, we devise a spatial regularization term f (Z) to incorporate the neighbor relationship information in X. The neighbors of the i-th column x i (e.g., x i−1 , x i+1 ) could be similar to x i . In other words, in the coding space, the corresponding neighbors of z i , which is a new representation for x i , could be close to z i . The spatial regularization term f (Z) is defined as: 2,1 -norm adopted here aims to achieve the similarity between neighboring pixels. W is the weight matrix enforcing the neighboring correlations in X, which is defined as: where ω is the spatial adjacent matrix of X, ω is the lower To show the spatial regularization term more clearly, denoted the adjacent relationship as ω, we illustrate a simple case in Fig. 3. Specifically, W is induced as follows: , which is devised to enforce the neighbors more probably to be segmented into the same superpixel. Locality-Constrained Subspace Clustering Model: Integrated with the spatial regularization in (2), we propose a locality-constrained subspace model for image superpixel segmentation: where λ 1 and λ 2 are the trade-off parameters to balance the weights of global regularization and spatial regularization, respectively.
Algorithm 1 Locality-Constrained Subspace Clustering Superpixel Input: An image I, desired number of superpixels K, initial number of pixel units n Output: K superpixels 1: Initially segment the image I into n raw pixel units by K-means. 2: Extract features of the raw pixel units as input data set X = [x 1 , x 2 , · · · , x n ] (each column is the feature vector of the corresponding pixel unit). 3: Construct matrices ω, D and W . 4: Solve the coding matrix Z according to Algorithm 2. 5: Build the affinity graph G with respect to Z according to (5). 6: Segment G with Ncut to generate K clusters. 7: Merge isolated pixels with the nearest neighbor clusters and update K.
After working out the coding matrix Z (See the detailed optimization algorithm in Section III-C), an affinity graph G can be constructed for segmenting the coding matrix Z and determining the subspace clusters. Similar to [34], the affinity graph G is defined as: This symmetrization transformation aims to assure the two connected points are in each other's sparse representation, to achieve the ideal correspondence of each element to the points from the same subspace. Then, Ncut [39] is adopted to generate the superpixel clusters with G. Moreover, in image superpixel segmentation, there may exist some isolated pixels not belonging to any clusters after the clustering step. Generally, such pixels are enforced to be assigned to the nearest cluster. The complete procedure of the proposed method is summarized in Algorithm 1. We show the result using our locality-constrained subspace clustering method in Fig. 2(f). It is clear that compared with other methods, only our method preserves the boundaries of the spots. More experimental results are presented in Section IV.

C. Optimization
To solve the objective function in (4), we adopt alternating direction method of multipliers (ADMM) [45] to devise an optimization algorithm. Two auxiliary variables U = Z and V = U W are induced to separate the terms of variable Z. Hence, (4) is equivalent to be written as: The augmented Lagrangian function of (6) is: where Θ ∈ R n×n and Ξ ∈ R n×n are the Lagrangian multiplier matrices, α 1 and α 2 are the regularization parameters. We can solve for U , V and Z alternatively by minimizing L when fixing the other two.
1) Update Z when fixing U and V: The problem in (7) becomes Since (8) is element-wise to each Z ij , by soft thresholding operator [46], we can get a closed-form solution: 2) Update U when fixing Z and V: By setting the derivative with respect to U to zero, we can obtain the following equation: (11) is a standard Sylvester equation [47], which can be linearized and solved by: where ⊗ is the tensor product and vec(·) returns a column vector by stacking the columns of a matrix.

D. Discussion
In the proposed model, three blocks of variables Z, V , and U need to be optimized in eq. (6). In each iteration, with the help of the ADMM, the objective function can converge to a local optimum with the optimized variables. Though there is no general convergence proof for the ADMM based algorithm with more than two blocks of variables [49], to the best of our knowledge, we observe that our algorithm empirically converges well and produces satisfactory performance on realworld data sets, contributing to the fact that the proposed model is convex, and each subproblem has a closed form solution.
The convergence curve is illustrated in Fig. 4. The algorithm can be updated iteratively until it the stopping conditions are satisfied, i.e., the iteration error is reduced to a stable value (10 −3 in our experiments) after several iterations or the algorithm exceeds the maximum iteration. We empirically find that 20 iterations are sufficient for most images, and thus we set the maximum iteration number to 20 for efficiency.
The computational complexity of the proposed method is dominated by the initialization by K-means and the optimization process. The computational complexity of K-means is O(N ), where N is the total number of pixels in the input image. That of the optimization process is O(k 3 ), where k is the initial number of pixel units, which is far less than the total pixel number N. Hence, the time complexity for the proposed method is O (N + k 3 ).

IV. EXPERIMENTAL RESULTS AND DISCUSSIONS
To show the merits qualitatively and quantitatively, we compare the proposed method with some state-of-the-art methods, including Ncut [39], Turbopixel [40], TPS [2], LRW [27], SLIC [28], ASS [31], FLIC [30], LSC [29] and BASS [41]. The parameters of the compared methods are set according to the original works. Since BASS cannot set the superpixel number manually, we adopt the published results directly.

A. Evaluation Criteria and Datasets
In the experiments, three commonly used metrics are utilized to evaluate all of the methods quantitatively: achievable segmentation accuracy (ASA), boundary recall (BR) and under-segmentation error (USE) [30]. ASA measures the highest achievable accuracy of a segmentation, i.e., the percentage of the most labeled pixels of superpixel segments overlapped from the ground truth segments with respect to the total number of the pixels in superpixel segments. Higher ASA indicates better segmentation of image content. USE measures the fraction of pixels of superpixel segments that not overlapped from the ground truth segments. Lower USE means fewer pixels leak from the object boundaries and edges. BR is the fraction of correctly overlapped pixels of superpixel segments from the ground truth boundaries. A high BR means the image content boundaries are well preserved after segmentation.
We compare the proposed method with some state-of-theart superpixel segmentation methods on the public Berkeley Segmentation Dataset (BSDS500) [50] for thorough superpixel segmentation evaluation. BSDS500 contains 500 natural images including humans, animals, buildings, landscape and outdoor scenes. Since the dataset has accurate annotated segmentation boundaries, almost all of the superpixel methods adopt BSDS500 for evaluation. The dataset provides 5 ground truth for each image and we average segmentation results with respect to the different ground truth for evaluation. For further comparisons, we evaluate the methods on two other datasets: Densely Annotated VIdeo Segmentation (DAVIS) [51] and Cityscapes [52]. Since these datasets do not contain accurate boundaries, the ASA metric is adopted to evaluate the upper bound segmentation performance. Besides, we add NYUv2 [53] dataset which consists of indoor-scene images for further extensive comparison.

B. Parameter Selection
The image size of BSDS500 dataset is 481 × 321. If every pixel is drawn as a sample, the dimension of the input is 154401. To speed up the computation, we use K-means to produce raw pixel units to keep the dimension of input data at thousand level. Let r denote the ratio of the initial number of pixel units n and the desired number of superpixels K. Fig.  5 shows the influence of r with respect to different metrics. We can see that the segmentation quality keeps a relatively high level when r is around 3. Indeed, blindly increasing the number of raw input units cannot improve the performance obviously, but increase the computation amount. Eventually we set the parameter r to 3 in our experiments to trade off the computation speed and segmentation performance. We have tried different λ 1 and λ 2 , and find that excessively increasing the weight of λ 1 and λ 2 may reduce the performance. Finally, λ 1 and λ 2 in (4) are empirically set to 10 −6 and 10 −4 , respectively (see more details in Section IV-E). The features are extracted from each pixel units, including color feature (average gray and R, G, B values of pixels unit), spatial feature (average coordinate of pixels unit), edge feature (average x and y components of gradient, gradient magnitude, and gradient direction) and texture feature (extracted by [54]). These features are concatenated and vectorized to X.

C. Visual Comparison
We compare our method with some state-of-the-art superpixel segmentation methods. Fig. 6 illustrates some visual comparative results when the desired number of superpixels K is set to 300. As can be seen, the proposed method maintains more detailed and complete object boundaries. For example, in the first image, our method has kept more detailed boundary information: most of the small spots in the butterfly wing are captured by our method, but almost unrecognizable by the compared methods. For the second image, our method also catches the most complete and accurate boundaries, and moreover, the homogeneous regions of the zebra stripes have better coherence and clearer segmentation compared with others. In the third image, the boundaries of doggies are preserved more accurately than other methods. Ncut, LRW and TPS cannot catch most of the boundaries in the doggies' forehead and nose, while SLIC, FLIC and LSC can capture most boundaries but missing a little bit. Only the segments of our method can cover all the boundaries of foreheads, and the most boundaries of noses. The reason is that in our method, each superpixel is regarded as a subspace containing pixels with similar intrinsic attributes, such as similar texture, boundary, color and coherent spatial information. By this way, the pixels in the unit containing independent semantic information are more probably segmented into one superpixel, no matter how small the unit area is.
It is clear that other methods can also capture more image boundaries like ours if the desired superpixel number K increases heavily. In Fig. 7, we investigate the minimum desired superpixel number K for accurate spot segmentation in the butterfly image. The proposed method can segment the spot when K = 300, while that of LRW, SLIC, and LSC are 1600, 1800, and 1600, respectively. This shows that our method can achieve more details and semantic information with less superpixels compared with other methods.
Subsequently, the proposed method is compared with three recent deep learning based methods, and Fig. 8 illustrates some visual results. It can be seen that our method can preserve more detailed boundaries than the deep-learning based methods. In the first image, only SEAL and our method segment the logo of the airplane. In the second image, both the two athletes' number "1257" and "282" are mostly recognizable in the proposed method, while in FCN, SSN, and SEAL, the number "1257" are unrecognizable. So our method is comparable with deep learning based methods, and even has superiority in detail preservation. It is also worth pointing out that our method is unsupervised, which does not need any ground truth information, while the compared deep-learning based methods all need some labeled samples that are costly to acquire.  Fig. 9 shows the quantitative comparisons with some stateof-the-art methods. Fig. 9 (a) and Fig. 9 (b) are the BR and USE metric results, which show that our method can largely improve the performance and outperforms all the other compared methods. Taking 300 superpixels for example, our method's minimum percentage gain (computed with the highest score of the compared methods) of USE is 62.5%, while that of BR is 6.7%. Fig. 9 (c), Fig. 9 (d), and Fig.  9 (e) are the ASA scores on BSDS500, DAVIS, Cityscapes datasets, respectively. In terms of ASA, the proposed method is competitive and improves the performance especially when the number of superpixels is small. Fig. 9 (g), and Fig. 9 (h) are the BR and USE score on NYUv2 dataset. The performance of our method also achieves the best performance in this dataset. The minimum and maximum percentage gains of BR are 3.60% and 41.15%, respectively. The minimum and maximum percentage gains of USE are 3.43% and 40.22%, respectively. The proposed method can largely improve the segmentation quality. Fig. 9 (i) shows the BR scores of our method and deep learning based methods FCN, SEAL and SSN on the 200 images of test in BSDS500. Compared with the deep learning based methods, the BR score of our method is higher than SEAL and FCN, and comparable with SSN. Moreover, our method performs better than SSN when the number of superpixels is less than 300.

D. Quantitative Comparison
Our method is implemented in Matlab and runs on a PC with a 3.6 GHZ Intel Core i7-6850K CPU and 16GB RAM. Fig.  9 (f) shows the trade-off between efficiency and performance compared with some state-of-the-art methods when desired superpixel number is 100 and 600. According to the BR    Fig. 10. Ablation studies on BSDS500. When λ 1 = 10 −6 , λ 2 = 10 −4 , BR, USE, and ASA can achieve the highest score.
scores, the proposed method can achieve the best segmentation performance, slightly slower than SLIC and LSC, but faster than others.
In conclusion, the quantitative results show that the pro-posed method performs better than the compared state-of-theart methods by standard evaluation criteria.

E. Ablation Study
To verify the effectiveness of the global regularization (GR) and spatial regularization (SR) in the proposed model, we evaluate our method on BSDS500 by varying the trade-off parameters λ 1 (the weight of global term in (4)) and λ 2 (the weight of spatial term in (4)). Table I lists the results of our model without the GR or SR respectively with different K. E is the error term in (4). It can be seen both the SR and GR are important to the proposed model, since E + GR + SR always performs the best. Moreover, E + SR performs better than E + GR, indicating the importance of the spatial prior term.
Moreover, Fig. 10 illustrates the BR, USE, and ASA scores with different λ 1 and λ 2 . From Fig. 10, we can see that excessively increasing the weight of λ 2 reduces the performance of BR. When λ 1 = 10 −6 and λ 2 = 10 −4 , all the metrics achieve the highest score. Therefore, we choose these two values for the parameters. For DAVIS and Cityscapes datasets, we adopted the same parameter values and achieved excellent results.

V. CONCLUSION
In this paper, we formulated image superpixel segmentation as a subspace clustering problem based on the assumption that an input image is constituted of several subspaces, and pixels with similar intrinsic attributes lie in the same subspace. We designed a novel local regularization term to enforce the spatial correlation, and the proposed spatially constrained subspace clustering model is able to generate content-aware superpixels with more detailed boundaries. The proposed method has been compared with some state-of-the-art methods, and the experiments on different standard datasets demonstrate that the proposed method achieves superior performance both quantitatively and qualitatively.


Hydrotreating unit models based on statistical and fuzzy information

This paper presents the results of mathematical models’ development for hydrotreating reactor, stripping соlumn, absorbers and hydrotreating furnace which are basic units of hydrotreating block in catalytic reforming unit. Since these objects of modeling in reforming unit of Atyrau refinery operate in conditions of insufficiency and fuzziness of the initial information, their mathematical models are developed on the basis of a systematic approach, using available information of different nature (experimental-statistical data, fuzzy information from the experts), with appropriate methods of construction for mathematical models. Mathematical models, describing the dependence of the production output from the hydrotreating reactor, columns and furnace, are designed as a nonlinear regression models based on experimentalstatistical data. Whereas models, evaluating the quality indicators of generated products from the hydrotreating reactor and columns, i.e., hydrogenation, hydrogen-containing and hydrocarbon-containing gases are built based on fuzzy information from specialists-experts in the form of fuzzy multiple regression equations. We have plotted the graph for the dependence of hydrogenation products output on the temperature in the hydrotreating reactor. To describe the dependence of the optimum temperature value in hydrotreating process on raw material quality, a linguistic model is designed based on compositional rules of inference and fuzzy information. Membership functions of fuzzy parameters are constructed for linguistic models.

Introduction
In this study, authors discuss in more detail the first direction of improvement for hydrotreating processes, based on hydrotreating unit of catalytic reforming installation in Atyrau refinery. There are known research works on mathematical modeling and optimization of technological objects and processes of oil refining, including hydrotreating process [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. However, in practice, production situations can often be associated with the of insufficiency and fuzziness of the initial information, therefore the problems of modeling and optimization of their operation modes, formulation and solutions using traditional methods do not provide adequate solutions These objects include hydrotreating unit LH in catalytic reforming of Atyrau refinery, basic units of which operate in the conditions of uncertainty associated with the randomness and fuzziness of initial information [19][20][21][22][23][24][25][26]. In order to create a mode optimization and control system of the hydrotreating unit in catalytic reforming facility, it is necessary to develop a system of mathematical models for main interconnected units of the facility, namely hydrotreating reactor R-1, steam column C-1, absorbers C-2 and C-3, and hydrotreating furnace F-101. In the development of mathematical models of these basic units of hydrotreating facility there may be a problem related to insufficiency, uncertainty and fuzziness of the initial information. The uncertainty can be caused due to deficiency, randomness and fuzziness of the available information, which is necessary for the development of mathematical models of objects. In these cases, one has to use an appropriate approach, such as a hybrid method proposed in [27][28][29][30][31][32][33][34][35][36][37][38], which allows constructing mathematical models of objects based on the available information of various character.
Let be and -a set of available input parameters of the object, which is characterized by the probabilityand fuzziness -values are determined using instrumentation, but are characterized by randomness. Values are fuzzy and estimated by human, i.e., specialists-experts based on their knowledge, experience and intuition fuzzily, i.e., by words and phrases. Based on such initial information of the different nature, it is necessary to identify the structure and parameters of models for main units of hydrotreating facility. Consequently, it is required to use the methods of probability theory, experimental and statistical modeling techniques, as well as known methods modified and adapted for operation in the fuzzy environment, for example the method of forward selection in regression, and the method of least squares [39][40][41][42][43][44][45][46][47][48][49][50], as well as a hybrid method of building models [51,52].
Thus, based on the aforementioned methods, it is necessary to develop a system of mathematical models for hydrotreating reactor R-1, columns C-1, C-2 and C-3, and hydrotreating furnace F-101. For R-1 it is needed to identify equations that enable calculating the values of the output parameters of the reactor: hydrogenate volume; and evaluating the values of quality indicators of the product: -respectively, unsaturated hydrocarbons, sulfur and water-soluble acid and alkali in the composition of hydrogenate. For columns C-1, C-2 and C-3 it is necessary to identify models that enable determining the hydrogenation output from C-1, HCG from C-2 and hydrocarbon-containing gas from C-3, as well as evaluating the basic quality indicators of column products from input parameters of columns. Mathematical models of the hydrotreating furnace F-101 must determine the amount of feed with gas and flow temperature from the exit of furnace depending on values of F-101 input parameters.
The aim of this work is to develop a system of mathematical models for basic units of hydrotreating block in Atyrau refinery catalytic reforming based on the available information of different nature, which can be used in optimizing the process parameters and modes control of hydrotreating unit.

Material and methods
Hydrotreating unit of catalytic reforming installation is one of the main units of the installation [53][54][55][56][57][58][59][60]. In the hydrotreating unit, hydrotreating process of straight-run gasoline from crude oil distillation occurs. Accordingly, there is increase in oil quality by removing the sulfur from their composition and as well as other harmful compounds and impurities that degrade the performance characteristics of the process equipment and metal units [61][62][63][64][65][66][67][68][69][70][71][72][73][74][75]. Thus, the process of hydrotreating reduces corrosion of metal equipment and pollution of environment and atmosphere. Hence, an actual task of technology science and petroleum production is the study and improvement of hydrotreating processes in refinery based on research-based methods, such as methods of mathematical modeling and optimization.
Consider the description of the technological diagram for hydrotreating unit in catalytic reforming of Atyrau refinery (AR) shown in Figure 1. The catalytic reforming process is designed to produce high octane motor gasoline, feedstock for petrochemical synthesis (aromatic hydrocarbons) and for hydrogen-containing gas (HCG), which are used in hydrogenation processes in refining. The target product of the catalytic reforming LH-35-11 / 300-95 in Atyrau refinery are high-octane commercial gasoline and liquefied natural gas. The installation LH consists of 4 units: preliminary hydrotreating unit of feedstock, i.e., straight-run gasoline, produced at the unit of primary oil refining; catalytic reforming unit, which increases the octane number of the hydrotreated gasoline; deethanizer and catalyzate stabilization unit; purification unit of circulating hydrocarbon-containing gases by monoethanolamine (MEA) and regeneration of MEA [76][77][78][79][80][81][82][83]. Feedstock for hydrotreating unit LH in Atyrau refinery is straight-run gasoline, produced in primary refining unit AVT ET-3. Hydrotreating process proceeds in HCG environment and it relates to a catalytic process, designed for purification and removal of organic compounds of sulfur, oxygen and nitrogen, which are poison for the catalyst, from the feedstock and from straight-run gasoline. Feed from reservoir tank is supplied with HCG by A16 pump for mixing. The mixture of feed and HCG is supplied to serially connected 3 heat exchangers HE-1/1-3, here due to the counterflow of aerated feed from the reactor R-1 and HE-3 reboiler is heated to a temperature of 260 °C, then further supplied to the hydrotreating furnace F-101. From furnace F-101, the mixture of feed and gas with a temperature of 300-343 °C is fed to hydrotreating reactor R-1. Hydrotreating reaction proceeds in the reactor R-1 with use of S-12 catalyst, i.e., feed is preliminarily hydrorefined from sulfur, nitrogen and oxygen containing compounds. Heat of the mixture of unstable hydrogenate and circulating gas from the reactor output and heat of the gas reaction at a temperature 340-420 °C are used for heating the mixture of feed and gas firstly in heat exchanger HE-3 of stripping column SC, then in heat exchangers HE-1/1-3 [84][85][86][87][88][89][90]. Products in the form of gas, after cooling to 35 о С temperature in refrigerators RR-101 and RF-1 are supplied to the separator S-1. HCG is separated from the liquid in S-1 and is fed to the absorber C-2 for removal of hydrogen. Gas from the exit of the absorber C-2 after passing through the separator S-4 is divided into two streams: 1) circulating gas, after compression in the compressor is fed back to the system of feedstock hydrotreating; 2) excess HCG from the exit of unit, the liquid phase of separator S-1 passes through the heat exchanger HE-2, where it is heated to a temperature of 150 о С and fed to 7, 9, 23 tray of stripping column C-1. Hydrogen sulfide and water are stripped from C-1 column, which has 30 trays, from hydrogenate at temperatures up to 270 °C and pressures up to 15 atm, in addition to that light hydrocarbon are discharged through the top of column.
After stripping column C-1, overall composition of sulfur compounds in hydrogenate must not exceed 0.0005 wt%. Gases in the vapor state from the top of column C-1 are out at a temperature of 135 °C, they pass through refrigerators RR-101 and RF-1 and fed into the separator S-2 with a temperature of 35-40 °C. From the separator S-2, it is fed back to stripping column C-1 in liquid phase. Precipitated water in the separator S-2 is discharged into the sewers. Hydrocarbon gas from the separator S-2 enters the absorber C-3 for purifying from hydrogen sulfide. Hydrocarbon gas from the top of the absorber C-3 is fed into the fractionating absorber C-6 or the fuel tank of factory. Thus, during hydrotreating process, chemical transformation of substances occurs under the influence of hydrogen gas at high pressure and high temperature [91]. During the hydrotreating process, sulfur compounds are reduced in petroleum products, additional unsaturated hydrocarbons are saturated, resin composition and oxygenate-containing are lowered and hydrocracking of hydrocarbon molecules are done. Improvement of hydrotreating processes in refinery using modeling techniques allows the following [92]: conducting hydrotreating processes in the optimal mode; improving the quality indicators of generated products.

Results and discussion
To solve this problem, i.e., to develop a system of mathematical models for the main units of the hydrotreating facility, which determine the dependence of the output parameters of the unit (product, quality) on the input mode parameters using available information of various character. For example, to develop mathematical models of hydrotreating reactor that enable determining the amount of hydrogenate from the output of the hydrotreating reactor R-1, authors use experimental-statistical data, characterized by probability and expert information, having fuzzy nature [93; 94; 95-110]. For structural identification of reactor R-1 models, the idea of forward selection in regression, whereas identification of the model parameters is based on a modified method of least squares. Thus, the mathematical models for the reactor R-1 of hydrotreating unit are designed using statistical data and fuzzy information processed by methods of mathematical statistics and the theory of fuzzy sets. Fuzzy information is collected and formalized by the methods of expert assessment and fuzzy set theory. Based on processing of experimental-statistical data and expert information, as well as by a method for constructing the fuzzy model [111; 112-121], authors conduct structural identification of models that describe the quality of products from hydrotreating reactor R-1 in the form of the following fuzzy equations of multiple regression: , where: -unsaturated hydrocarbons in the composition of the product, i.e. hydrogenate (not more than, i.e. 1%); -sulfur in hydrogenate content ( 0.00005%); − water-soluble acids and alkali in hydrogenate content ; -feed, in this case, straight-run gasoline (45-80 m 3 /hr.); − reactor pressure (20-35 kg/cm 2 ); − reactor temperature (300-343 о С); − rate of feed supply (0.5-5 h -1 ); − circulating hydrogen-containing gases (HCG) -ratio hydrogen/ hydrocarbons (200-500 nm 3 ); , -= identified fuzzy coefficients, free term; linear effect; quadratic effect and mutual effect, respectively. The allowable values of the output fuzzy parameters, and intervals of changing input and operating parameters are denoted in the parenthesis [122]. To identify the unknown parameters (regression coefficients) of the model (1): and -membership functions of fuzzy sets describing the quality of hydrogenate are distributed to the following level set α: α=0.5; 0.85; 1. Since in this case the membership function has the bell-shaped form, values of fuzzy parameters are obtained at 5 points α=0.5; 0.85 (left); 1; 0.85; 0.5 (right). Values of the input, regime and output parameters are observed for each selected α level. Thus, the obtained models, describing the product quality from the exit of reactor R-1 in the form of multiple regression for each α level. Since the obtained equations have a form of regression equations, the problem of identification of the unknown coefficients can be resolved using known parametric identification methods, for example by the method of least squares. In this work to identify the regression coefficients authors use REGRESS program package, which is based on modified least squares method to determine the regression coefficients of the linear and nonlinear regression models with any number of input parameters . In this manner, after parametric identification, mathematical models, describing the effect of input, mode parameters on the quality of the hydrogenate, i.e., on the content of unsaturated hydrocarbons , sulfur and water-soluble acids and alkalis for each α level have the following form: The identified values of the coefficients are combined together with the help of the expression below from the fuzzy set's theory [122]: or (5) where: .
In the derived models, regressors which have no effect on . and or very little influence are set to zero, i.e., not shown. As a result of conducted study and data processing data it is determined that for finding the volume of product output from the reactor R-1, i.e., the amount of hydrogenate -on the basis of experimental-statistical data, one can construct a statistical model, which allows estimating the value from the input and mode parameters using non-linear regression equations. After identifying the structure and parameters of this model, in the same way to the above approach, a mathematical model that allows calculating the volume of hydrogenate from the exit of reactor R-1 has the following form: Authors have constructed the plot describing the dependence of hydrogenate output on temperature in the reactor for the fixed input values of feed input and the other mode parameters: and . (Figure 2). The influence of the other mode parameters and and their mutual influence on the output parameters of the reactor are considered. The linguistic model is built to determine the optimum temperature of hydrotreating process on the basis of expert information, compositional rules of logical inference and base of rules. The resulting linguistic model has a logic relationship: "If the feed is heavy, then the temperature of process is low, otherwise, if the feedstock is light, the temperature of process is high". On the basis of expert assessment and using methods of fuzzy sets theory in the form of an exponential dependence, the membership functions are constructed that describe the fuzzy parameters of linguistic model: heavy (low thermal resistance) feed -straight-run gasoline: ; light feed: ; (9) low temperature: ; (10) high temperature: .
Thus, the structure of linguistic model, evaluating the optimum temperature in dependence on quality of feed: , is determined by compositional rules of logical inference: .
In the obtained linguistic models, the following notations are introduced: ts, ls, nt, vt -"heavy feed"; "light feed"; "low temperature" and "high temperature", respectively; , -input and output linguistic variables, which describe the quality of feed and the optimum temperature of hydrotreating process, respectively; ,fuzzy subsets that describe and . Let us provide the results of the development of mathematical models of columns C-1, C-2 and C-3, which are also the main units of hydrotreating facility. Stripping column C-1 of hydrotreating unit is intended to separate water vapor and hydrogen sulfide from the product, i.e., of hydrogenate. Content of sulfur compounds in the product must not exceed 0,0005 wt%. Column C-2 of hydrotreating unit is absorber, which is intended for purification from hydrocarbon gases. Whereas, column C-3 of hydrotreating unit is also absorber, where purification of hydrogen sulfide from hydrocarbon gases occurs. Based on study results of operating modes of columns C-1, C-2, C-3, and data analysis, as well as taking into account additional fuzzy information provided by the specialists-experts to develop mathematical models of these columns, it is decided to use a hybrid method of modeling [123; 124]. Based on research work results of columns C-1, C-2 and C-3 and the results of expert assessment, authors identify the following input (x) and output (y) parameters describing the mode of operation: volume of feed at the inlet of the columns; inlet temperature; pressure in the columns; supply rate in column C-1; volume of hydrogenate from the exit of column C-1; volume of HCG from the exit of column C-2; volume of hydrocarbon-containing gases from the exit of column C-3; content of sulfur compounds in hydrogenate, quality of C-1 product; content of HCG, quality of C-2 product; content of hydrocarbon-containing gases, quality of C-3 product. Input parameters for columns , i.e., feed consumption, inlet temperature, pressure in the columns and the supply rate for column C-1 and the output parameters describing the volume of the product at the output of columns are measurable parameters. This means that it is possible to collect statistical data of these parameters: Whereas, quality indicators of products generated in columns K-1, K-2 and K-3 are described fuzzily. Therefore, the fuzzy quality indicators are estimated by experts, formalized and processed using methods of fuzzy set theory. Thus, the structure of mathematical models of stripping column C-1, absorbers C-2 and C-3 in hydrotreating facility is identified by the method of forward selection in regression and its modifications in the form of the following nonlinear regression and fuzzy regression equations: Identification of parameters, i.e., regression coefficients of regression models (14) is performed based on the method of least squares using statistical data of object. To identify the fuzzy parameters of fuzzy regression models (15) based on level set α, fuzzy equations are converted to a system of non-fuzzy equations, which are equivalent to the original fuzzy equations. Then similarly to the procedure of identification of fuzzy parameters for reactor R-1 model (1), fuzzy coefficients for different levels of α are identified through the method of least squares criterion to minimize the error between the model and experimental (real) data. To determine the volume of target product from column C-1, i.e., hydrogenate after the parametric identification of the model (14) for = j 1, the following results are obtained: (16) Also, with method of forward selection in regression and least squares, structures (14) and parameters are identified for absorber columns C-2 and C-3, i.e., volumes of HCG ( ) and hydrocarbon-containing gas ( ) are determined from absorber C-2 and from the column C-3, respectively: ; (17) . (18) Similarly, to the model (1) for the assessment of quality indicators of hydrogenate , it is possible to identify the fuzzy parameters of model (15) for evaluation of the product quality from columns K-1, K-2 and K-3: . Mathematical models for hydrotreating furnace F-101 of hydrotreating unit. The cylindrical hydrotreating furnace F-101 is intended for heating the hydrotreating product, i.e., hydrogenate to the required temperature. According to the study and analysis results, the following main parameters are identified that influence the operation of the furnace F-101 and hydrotreating process: flow rate, volume of feed at the inlet of the furnace F-101 in the range of 60-80 m3 / hour; temperature at the inlet of furnace F-101, in the range of 170-190 °C; pressure in furnace F-101, in the interval 40-43 kg/cm 2 . As a result of analysis for the collected data, and operation modes study of the hydrotreating furnace, the experimental-statistical method has been selected for the development of its model. Optimum operation mode of furnace can be selected on the basis of the mathematical model, which describes the effect of input variables on the output parameters, i.e., provides information on the thermal operation of the furnace. Mathematical description, which is the basis of the mathematical model, must determine the thermal operation parameters of furnace [125; 126]. The main disadvantage of hitherto used methods for calculating the furnace characteristics is that these methods only determine the integral indices of the heat exchange process, and do not regulate the possibility of the furnace tube heating. In the recent time, the researchers have proposed modeling methods based on theoretical studies to determine local heat transfer values, for example, the zone method. In mathematical terms the idea of the zone method calculation is the replacement of integral-differential equations describing the process of heat exchange with the algebraic equations approximating them with the limited system. By solving the obtained algebraic equations, energy characteristics of heat transfer is defined, i.e., temperature and flows of local zones. For this purpose, furnace research is divided into the limited number of area and volume zones, which have the same radiation properties. This approach to the calculation of the furnace characteristics can provide sufficient accuracy, to improve the accuracy it is necessary to increase the number of zones. However, this process is quite complex and the collection of the necessary information for its application in practice is also difficult.
Simple models are required to model the operation of industrial furnaces in interactive mode and to enable fast acquisition of necessary information and results. For this reason, the analytical method of N.I. Belokon can be used as a basis for modeling algorithm, this method is based on the simultaneous solution of heat transfer and heat balance equations [126]. Regression models are identified to calculate the output parameters of the hydrotreating furnace F-101 based on experimental and statistical data. The law of distribution of random measurements is taken close to the normal distribution, i.e., . Thus, the structure of the model, evaluating output of hydrotreating furnace: volume of mixture of feedstock and gas, and temperature of the outlet flow from the furnace, is identified by the following nonlinear regression equations: (19) The following notation in the model (19) can be adopted: model parameters, which will be identified, the known method of least squares can be used for their evaluation; − furnace F-101 mode parameters, volume of feed ( ); inlet temperature of the furnace ( ) And pressure in the furnace ( ), respectively. The following are the results of identification of regression coefficients in the model (19) using processed statistical data and using REGRESS program: As a result of system analysis and expert assessment, it has been determined that development of statistical and deterministic models using appropriate methods is the most effective for heat exchangers and separators in hydrotreating unit, respectively. The proposed approach to the development of mathematical models of interrelated technological objects on the basis of the available information allows the development of various models of real technological objects in the conditions of insufficiency and fuzziness of the initial information. The structure of developed mathematical models for main units of hydrotreating facility: reforming reactor R-1, columns C-1, C-2 and C-3, reforming furnace F-101 is identified as non-linear regression equations. The equations, describing the volume of production from units, have the form of conventional multiple regression equations. Whereas, the equations, describing the quality of products from the main units (content of unsaturated hydrocarbons − , sulfur − and water-soluble acids and alkalis − ; quality of products from columns C-1, C-2 and C-3: ) have the form of fuzzy multiple regression equations. In conditions of fuzziness in both input and output parameters, i.e., when the input and output of hydrotreating reactor are described using linguistic variables, it is proposed to construct the linguistic models based on compositional rules of logical inference. This approach is implemented in constructing the linguistic model describing the dependence of the optimum temperature of the hydrofining process on quality indicators of feedstock. In this case, for the construction of the membership function of fuzzy parameters, it is necessary to select the exponential dependence, which has adjusting coefficients to more adequately describe the function. The generated models can be used when optimizing the process parameters, selecting optimal operation modes and controlling hydrotreating process.

Conclusion
This paper provides the research results of hydrotreating process fundamentals, the primary directions of modernization and improvement in refinery hydrotreating are highlighted. The main results are demonstrated in the direction of the hydrotreating processes in the optimal mode based on design and modelling of operation modes of hydrotreating reactor. The mathematical models of the main units in hydrotreating facility of reactor R-1 are obtained; stripping column and absorbers C-1, C-2, C-3; reforming furnace F-101 of catalytic cracking facility in Atyrau refinery are characterized by insufficiency and fuzziness of the original information. To solve the problems with fuzziness of initial information and the development of mathematical models, it is proposed to use the available information of different nature using the hybrid method of model development. The mathematical models of the main units in hydrotreating facility are based on experimentalstatistical data and fuzzy information from specialists-experts. The mathematical models for determining the volume of product from the exit of units are identified as statistical models of regression type. Whereas models, evaluating fuzzily described quality indicators of generated product, are defined in the form of fuzzy equations. Structural identification of developed models is conducted on the basis of the idea of forward selection in regression, whereas parametric identification is implemented by the modified method of least squares using REGRESS software package. The graph is plotted for dependence of the hydrogenate output on temperature in the hydrotreating reactor R-1 for fixed values of the other mode parameters. Using linguistic compositional rules of inference rules, the linguistic model is built that describes the dependence of the optimum temperature value on the feedstock quality.