Subpixel Mapping for Remote Sensing Imagery Based on Spatial Adaptive Attraction Model and Conditional Random Fields

For subpixel mapping (SPM), ensuring the operational efficiency of the algorithm and mitigating the effect of abundance errors often cannot be achieved simultaneously. To solve the problem, we propose a new SPM method based on the spatial adaptive attraction model (SAAM) and conditional random fields (CRFs). First, the proposed SAAM obtains the spatial adaptive attraction value by adaptively adjusting the spatial attraction value obtained using the traditional spatial attraction model, thereby turning the display form of the abundance constraints in the SPM into an implicit form for expression, to perform the physical significance of the abundance constraints with the relative size of the attraction value of each subpixel. Second, the spatial adaptive attraction value of the implicitly represented abundance constraints and the local spatial smoothing prior are modeled in the CRFs, and the model makes full use of the spatial information in the label field while considering the abundance constraint. Third, Graph-cut is used to optimize the model, the proposed SPM can not only guarantee the operational efficiency, but also extinguish the influence of abundance error and decrease the noise artifact on the results of SPM. Experiments on three remote sensing images show that the proposed SPM accuracy is considerably better than the previously available SPM methods and is the least time-consuming. This study provides a new solution for the SPM of remote-sensing images.


I. INTRODUCTION
L AND cover information obtained by classifying remote sensing images plays a vital role in various applications, such as environmental disaster damage assessment, precision agriculture, and urban planning [1], [2], [3]. Traditional hard classification methods treat each pixel in the remote sensing image as a unit by assigning class labels to each pixel while the spatial resolution remains constant. However, using traditional hard classification methods to extract land cover information accurately is difficult because of the mixed pixels containing multiple land cover classes. Mixed pixel is a common phenomenon in remote sensing images [4], [5], particularly in hyperspectral images. Spectral unmixing or soft classification techniques are widely used to solve the problem of mixed pixels [6]. However, the soft classification technique only estimates the abundance map and the area proportions of different land cover classes within the mixed pixel. It does not provide information on the spatial distribution of land cover classes within the mixed pixel.
The concept of subpixel mapping (SPM) was proposed by Atkinson in 1997 to obtain the spatial distribution information of various land cover classes within the mixed pixel in the remote sensing images and solve the land cover map of the subpixel scale [7]. SPM divides low-resolution pixels into high-resolution subpixels and assigns class labels to these subpixels according to certain criteria. SPM techniques, often considered the further processing of soft classification, transform the abundance maps obtained from soft classification into land cover maps on the subpixel scale. Therefore, SPM divides low spatial resolution pixels into several highly refined subpixels to produce a hard classification land cover map with high spatial resolution, which is highly valued in coastline mapping, land cover change monitoring, urban building recognition, and other applications [8], [9], [10], [11], [12], [13], [14].
Different SPM methods have been developed since the concept of SPM was proposed. As shown in Fig. 1, these SPM methods can be roughly divided into two types: the two-step strategy method and the optimization strategy method [15]. The two-step strategy SPM method consists of two steps. subpixel soft attribute values for all classes of all subpixels are estimated with the optimization strategy. Second, the class allocation for the subpixels is performed according to the subpixel soft attribute values and abundance constraints. The two-step strategy SPM method strictly adheres to abundance constraints, where the number of subpixels of each land cover class within the mixed pixel corresponds to the input abundance map. Examples of the two-step strategy SPM method include traditional subpixel/pixel spatial attraction model (SASPM) [16], pixel-swapping algorithm (PSSPM) [17], linear optimization model [18], and the spatial interpolation algorithm [19], [20]. In addition, some AI-based methods exist, such as genetic algorithm [21], [22], BP neural network algorithm [23], and artificial immune algorithms [24]. Abundance constraints are valuable and meaningful in determining the spatial distribution of subpixels for different land cover classes, as the abundance map contains some useful information about the land cover classes. The less the abundance map error, the more useful the information. However, soft classification or spectral unmixing remains an open problem. The obtained abundance map does not necessarily accurately predict the true proportion of land cover classes within the mixed pixel, thereby affecting the performance of SPM techniques.
The optimization strategy SPM method was developed to reduce the effect of abundance error caused by soft classification or spectral unmixing on SPM. This type of SPM method relaxes the constraints on abundance and eliminates the uncertainty caused by soft classification. Furthermore, the optimization strategy SPM method takes the SPM problem as an optimization problem, iteratively generating the optimal SPM results using different objective functions and optimization strategies. Examples of optimization strategy SPM method include the SPM model based on Markov random field [25], [26], the Hopfield neural network SPM model [27], [28], SPM model based on maximum a posteriori (MAP) [29], and spectral-spatial constraint model [30]. These methods have achieved good results in alleviating the effect of abundance errors on the result of SPM.
Although the optimization strategy SPM method can handle the impact of abundance error on SPM to a certain extent, its operation speed is often slow. Moreover, the optimization is very time-consuming and cannot handle remote sensing images with large regions. By contrast, although the two-step strategy SPM method is affected by the abundance error, it runs fast and can handle an extensive range of remote sensing images. The SPM method based on the spatial attraction model has obtained extensive engineering applications because of its high efficiency and clear physical relevance of spatial correlation. This study aims to solve the many noise artifacts in the result of the spatial attraction model due to the abundance error while maintaining high efficiency. Therefore, a new SPM for remote sensing images based on spatial adaptive attraction model (SAAM) and conditional random fields (CRFs) is proposed in this study with the following contributions.
1) An SAAM is designed, which obtain spatial adaptive attraction value by adaptively adjusting the spatial attraction value obtained using the traditional spatial attraction model. Spatial adaptive attraction values turn the display form of the abundance constraints into an implicit form for expression to maintain the relative size of the spatial correlation of each land cover class and characterize the physical significance of the abundance constraints. Thus, the smoothing of the noise artifacts is facilitated. 2) Spatial adaptive attraction value and the local spatial smoothing prior are modeled in the CRFs. The unary potential term designed from the spatial adaptive attraction values can calculate the cost that each subpixel is given the corresponding class label. The pairwise potential term modeled using a multilevel logistic (MLL) model favors the adjacent subpixel to adopt the same class labels, so that the spatial information in the label field can be fully utilized while considering the abundance constraint.
3) The CRFs model is optimized using the Graph-cut, make the proposed SPM can not only ensure the advantages of fast running speed and high efficiency of the two-step strategy SPM method, but also extinguish the influence of abundance error to a certain extent and decrease the noise artifact on the results of SPM. The rest of this article is structured as follows. The traditional spatial attraction model is described in Section II. The proposed SPM is described in Section III. Section IV conducts the experimental analysis on four remote sensing images. Section V has the relevant discussion. Section VI comprehensively concludes this study.

II. SPATIAL ATTRACTION MODEL
The SPM method based on the spatial attraction model was proposed by Mertens in 2006 [16]. This method quantifies the spatial correlation as a subpixel within the central mixed pixel subject to the spatial attraction of different land cover classes in the neighborhood mixed pixel. The spatial attraction values of each subpixel within the central mixed pixel are calculated from the abundance of each land cover class of its neighborhood pixels and the Euclidean distance between the subpixel and the adjacent pixels. After the spatial attraction value of the subpixel is obtained, class labels are assigned to the subpixels based on the abundance of each land cover class within the central mixed pixel. The neighborhood pixel of the central mixed pixel contains at most eight neighborhood systems, namely, the 3 × 3 neighborhood system of the central mixed pixel. The remaining pixels cannot be attractive because of the long Euclidean distance from the central pixel.
Suppose that the set of pixel sites for the whole remote sensing image is I, and the number of pixels is N . In this case, the pixel at site i in the image is represented as x i . The set of subpixel sites within each pixel is represented by J, where each pixel contains d 2 subpixels, with d being the scale factor of the SPM. The class label of the subpixel at the site j is represented as l i,j within the ith pixel. Therefore, the remote sensing image can be represented as X = {x i |i ∈ I}, and the labels of subpixels can be represented as L = {l i,j |i ∈ I, j ∈ J}. Moreover, the total number of land cover classes included in the image with K and each land cover class is represented as k(k = 1, 2, . . . , K). The spatial attraction values can be expressed by the following equation: where at k i,j represents the spatial attraction value of the subpixel at site i, j that belongs to class k, and the set is expressed as AT = {at k i,j |i ∈ I, j ∈ J, k ∈ K}. N p is the set of the neighborhood pixel sites for the pixel x i , and f k t is the abundance where the neighborhood pixel x t belongs to the class k. The abundance set is expressed as F = {f k i |i ∈ I, k ∈ K}. D t i,j is the Euclidean distance between the subpixel at site i, j and the neighborhood pixel at the site t. The schematic diagram is as follows. Fig. 2 defines the 3 × 3 neighborhood system and the coordinate system of subpixel according to the scale factor d = 3, where the upper left corner with the coordinate of (0, 0) is the starting point. Given that the center of the pixel and subpixel is the calculation point, the coordinate of the first pixel in the upper left corner is (1.5, 1.5), and the coordinate of the first subpixel is (0.5, 0.5). Assuming that the coordinate of the pixel is (P a , P b ) and the coordinate of the subpixel is (S a , S b ), the distance between the pixel and the subpixel is calculated as follows: The spatial attraction value is the soft attribute value of the subpixel belonging to a certain class. The larger the value, the greater the possibility is that the subpixel belongs to this class. After the spatial attraction values are obtained for all the classes of each subpixel within the central mixed pixel, the continuous soft attribute values are converted into discrete subpixel class labels according to the abundance constraints and class allocation process Equation (4) satisfies two of the class allocation requirements. The first one is that, for each subpixel, it should be assigned to one and only one class. The second is that for each class, the number of subpixels belonging to it should be consistent with the abundance, and they should be completely exhausted during the class assignment, namely the abundance constraints. The abundance constraint specifies the number of subpixels belonging to the different land cover classes within the central pixel, thereby presenting the physical significance of the SPM based on the linear mixed model. The abundance constraint for (3) is display form, which enforces the number of subpixel of each land cover class in the mixed pixel according to the abundance and the scale factor. If the label of the subpixel is directly obtained according to the class corresponding to the maximum spatial attraction value of the subpixel, the proportion of each class within the mixed pixel obtained by spectral unmixing cannot be satisfied However, if the subpixel label is assigned strictly according to the abundance constraint, the abundance error can have a considerable impact on the performance of the SPM, namely, the noise artifacts in the SPM results that caused by the abundance error.

A. Spatial Adaptive Attraction Model (SAAM)
To satisfy the abundance constraint is an important problem in SPM, and it is the key to maintain the significance of the SPM algorithm. The abundance constraints of traditional spatial attraction model is in display form, which limits its further combination with other models (e.g., CRFs). In order to facilitate the integration with the CRFs, thus allowing for the spatial smoothing of the noise artifacts, and weakening the influence of the abundance error on the results of SPM, the SAAM is proposed in this part by referring to the adaptive visiting order of classes method proposed in [31] based on the traditional spatial attraction model. The objective function is shown as follows: where aat k i,j represents the spatial adaptive attraction value of the subpixel at site i, j that belongs to class k, and the set is expressed as AAT = {aat k i,j |i ∈ I, j ∈ J, k ∈ K}. The model adjusts the display form of the abundance constraints of spatial attraction value to an implicit form with an adaptive scheme, while maintaining the relative size of the spatial correlation of each land cover class. That is, when the label of the subpixel is directly obtained according to the class corresponding to the maximum spatial adaptive attraction value aat k i,j of the subpixel, the proportion of each class within the mixed pixel obtained by spectral unmixing can be satisfied, as shown in (7). So, the spatial adaptive attraction value satisfies the requirements of integration with CRFs. The calculation steps are as follows: Step-1: The pixel x i is selected and the number of subpixels for each land cover class within this pixel is calculated , k ∈ K based on the abundance f k i , k ∈ K and the scale factor d. According to (10), the local Moran index R k i , k ∈ K of each class of the pixel x i is calculated. Meanwhile, the d 2 subpixels within the pixel x i are marked as 0 value, namely Step-2: The local Moran index R k i , k ∈ K of the pixel x i is sorted by the value, and the visiting order of the class of the pixel x i is determined by that order, that is, the class k with the large local Moran index, whose corresponding at k i,j , j ∈ d 2 will be visited first.
and can be expressed as at The is the sequential priority, independent of the size of the at k i,j value, and (1), (2), . . ., (K) represents the order of class visited). Step Assuming the current visited class k = 2, the scale factor d = 3, and the C 2 i = 3, then the result of the descending order of the spatial attraction values belonging to the class 2 of x i can be expressed as 6 , and the first three subpixels are marked as

Algorithm 1: SASM.
Input: the abundance F , the spatial attraction value AT , the number of pixels N , the number of class K, the scale factor d, and the local Moran index R k i , k ∈ K, i ∈ N Output: the spatial adaptive attraction value AAT Rank the at

for end for end for
Step-4: All the spatial attraction value at k i,j , j ∈ d 2 of the current visited class k are adjusted according to (8) to obtain the spatial adaptive attraction value aat k i,j , j ∈ d 2 : (8) Following the class visiting order determined by step 2, the d 2 spatial attraction value at k i,j , j ∈ d 2 of the class k for the next visit of the pixel x i are treated according to (9). The adjusted spatial attraction value at k i,j , j ∈ d 2 are then recycled to Step 3 Step-5: All of the pixels in the image do this operation. Finally, spatial adaptive attraction value aat k i,j is output. The complete process is summarized in Algorithm 1.
The proposed SAAM has the following features. 1) Compared with the display form of the abundance constraint of (3), the spatial adaptive attraction value adjusted by (8) achieves the effect of an implicit form of abundance constraint. For subpixels marked as 0, where M i,j = 0, (8) automatically suppresses the spatial attraction value of the currently visited class k to less than 0, and the spatial attraction value of the next visit class k remains unchanged. For subpixels marked as 1, where M i,j = 1, (9) automatically suppresses the spatial attraction value of the next visit class k to less than 0, and the spatial attraction value of the currently visited class k remains unchanged. This ensures that the number C k i of subpixels for the next visit class k is less than the number of subpixels marked as 0, and then, when the spatial attraction value of class k is operated in descending order, the M i,j corresponding to the first C k i spatial attraction values are all equal to 0. And so that only one of K spatial adaptive attraction value for each subpixel within the mixed pixel is the maximum and is positive, and the rest of the spatial adaptive attraction values are less than 0. Moreover, the number of subpixels of class corresponding to the maximum spatial adaptive attraction value of all subpixels in the mixed pixel satisfies the abundance constraint, namely (7), as shown in Fig. 3(a). According to the above, it can be concluded that, when the label of the subpixel is obtained according to the class corresponding to the maximum spatial adaptive attraction value of the subpixel, the two requirements of the class allocation are satisfied, namely (4). Therefore, the spatial adaptive attraction value implicitly expresses the abundance constraint, and through the adjustment of (8), the spatial adaptive attraction value of each subpixel in the mixed pixel can not only maintain the relative size of the spatial correlation of each land cover class and characterize the physical significance of the abundance constraints. Meanwhile, the adjustment in (8) and (9) is an adaptive and simple scheme that does not require manual participation.
2) The class visit order in the SAAM process affects the results. Fig. 3(b) shows the class visit order of two cases, which produces different results. Therefore, the class visit order is solved by the local Moran index to determine it reasonably. The local Moran index is based on the pixel, and the order of visits to each class within the pixel is determined based on the spatial autocorrelation of the local region of the pixel. As the center, the pixel defines a window of N w × N w size. In this study, N w takes the value of 3. Then, the spatial autocorrelation of the pixel (e.g., x i ) for the class k in the local region is quantified by the Moran index, as shown in the following equation: where f i n 1 and f k n 2 represent the abundance of the class k of any pixel within the local window centred on x i . f k i represents the mean of the abundance of all class k within the local window.
The class with large R k i values is prioritized by comparing R k i to specify the class visit order within x i . The local Moran index determines the order of class visits by considering the spatial correlation within the class without any a priori information. Thus, the process of the solution is adaptive.

B. Modeling of SPM Based on CRFs
In the Bayesian framework, SPM is usually a posterior probability distribution problem. The posterior probability of subpixel class label given the abundance of each class can be equally written in the following form by using the Bayesian rule [32]: (2) .
For terms (1) and (2) in (12), two special probabilistic frameworks can be used to model the posterior probability distribution: discriminative CRFs framework and generative Markov random fields (MRFs) framework [33].
The MRFs framework is a generative model that models the joint probability distribution of random variables. However, in an image processing task, the random variable is modeled by the posterior probability distribution and derived as a combined form of the prior probability distribution and the likelihood probability distribution, expressed as follows: where the label field L has the MRFs property, and the prior probability distribution p(L) can be modeled using the Gibbs distribution. p(F |L) is the likelihood probability distribution for the observed data. If a sufficient correct likelihood function can be designed to make the most accurate expression of the generation process of the observed data, then the MRFs framework can fully mine the information of the observed data theoretically. For the feasibility and convenience of computation and modeling, the MRFs framework assumes that the observed data are conditionally independent and that the observed data of individual or local regions are modeled. Thus, the likelihood probability distribution can be expressed in the form of a factor product. The prior probability distribution is modeled using MRFs in the generative MRFs framework. Thus, the spatial correlation information can be exploited. However, the likelihood probability distribution prevents the generative MRFs framework from exploiting the spatial correlation information of the observed data because of the independence assumption. Therefore, on the basis that the label field L has a MRFs property, Professor Lafferty of Carnegie Mellon University defines the discriminant CRFs framework to model directly the posterior probability distribution as a Gibbs distribution with the following form [34]: where Z is the normalization constant, C is a clique, and V C (L) is the potential function of the clique C. Equation (14) has no strict independence assumption. Thus, the requirement of several independent hypothesis distributions of the observed data in the MRFs framework can be overcome, and the spatial correlation information of the observed data can be modeled flexibly. The CRFs framework does not need to describe and express the observed data, avoiding the explicit modeling of the likelihood probability distribution in the MRFs framework. It is also directly focused on L itself, which simplifies the problem to some extent. The CRFs framework is a discriminant model that directly obtains the posterior probability distribution of the label field L in the SPM. According to (14), the corresponding energy function can be defined in the following form: According to the Bayesian MAP, the SPM label field L can be estimated by the maximum posterior probability: The SPM problem is converted by (16) to minimize the energy function E(L). In this study, only unary potential and pairwise potential are considered. The energy function E(L) can be expressed in the following form: where Φ i,j is the unary potential term, N i,j is the site set of the neighborhood system of l i,j , and Φ i,j i,t is the pairwise potential term computed over the neighborhood system of l i,j . λ is the tuning parameter of the unary potential and pairwise potential.
Unary Potential: The unary potential term is the modeling of the cost penalty for the subpixel to obtain a certain class label. It is related to the probability that subpixel belongs to each class; that is, the larger the probability is that the subpixel belongs to a certain land cover class, the smaller the cost penalty of assigning the subpixel is to the class label. This study uses the proposed spatial adaptive attraction value to model the unary potential energy in the following form: where p(l i,j = k) indicates the probability that the label of the subpixel at site i, j belongs to the class k. AAT is the spatial adaptive attraction value described above, derived from the abundance F of each class. Based on the SAAM, p(l i,j ) maintains the magnitude of the spatial correlation of the subpixel and implicitly characterizes the limits of abundance constraints.
Pairwise Potential: The pairwise potential energy is the interpretation of the spatial prior knowledge of land cover classes, where adjacent subpixel tends to have the same class labels that have important information in determining the land cover classes of the subpixel. This study models the pairwise potential energy function of the 8-neighborhood system using an MLL regression model with the following formula: where η measures the weight of labeling smoothing potential. As for the parameter η, we simply set it to 1 [35]. The pairwise potential energy based on (19) encourages adjacent subpixels to acquire the same class labels and punish nonadjacent subpixels for having different class labels, thereby achieving a smooth effect on the result.

C. Optimization of SPM Based on CRFs
In this article, we use Graph-cut based α-expansion inference algorithm [36][37] to optimize and solve the established energy function [i.e., (17)]. Graph-cut based α-expansion algorithm is an optimization algorithm for solving the energy function minimization problem, which has been verified to be a fast and effective algorithm in practice and experiments, so, it has received wide attention in various applications [38].
Graph-cut algorithms are usually applied to binary labeling problem and can quickly converge to a global energy minimum. However, SPM is usually a multiclass labeling problem. For the multiclass labeling problem, the Graph-cut based α-expansion designs a special local search algorithm that can be used for the energy function minimization problem of multivalued variables. The local search algorithm repeatedly calculates the global minimum of the binary labeling problem in its inner loop using the Graph-cut algorithm. That is, at each α-expansion step, the binary labeling operation is performed so that each subpixel is either kept the current label or converted to a specific label α ∈ K. This operation for subpixel labels is performed simultaneously, and thus is exponential number for any specific label of possible transformation, which ensures that the algorithm has a strong local minimal property. Therefore, the Graph-cut based α-expansion can reduce the multiclass labeling problem to a sequence of optimization subproblems with binary variables, which can be easily optimized by Graph-cut methods.

D. Complete Algorithm
The newly developed SPM (SACRF) framework based on SAAM and CRFs is summarized in Fig. 4.

IV. EXPERIMENT AND ANALYSIS
Four remote sensing imagery were selected in verifying and analyzing the proposed algorithm to verify the performance of the proposed method. The four remote sensing images include a simulated hyperspectral imagery, a synthetic hyperspectral imagery, and two real multispectral imagery.
The methods for comparison include the traditional PSSPM, SASPM, spatial-spectral interpolation algorithm (SSSPM) [39], and radial basis function interpolation algorithm (RBF-SPM) [40]. All the abundances involved in the experiment were obtained by performing a nonnegative least squares solution of the original images. A hyperparameter search with interval 1 was performed in the range of 1−20 to determine the value of the optimal tuning parameter λ. The results were quantified using the Kappa coefficient, overall accuracy (OA), and producer accuracy (Prod.Acc).

A. Experiment 1: Simulated Hyperspectral Imagery
Hyperspectral imagery was simulated using a linear mixing model [41] and zero-mean i.i.d. Gaussian noise [42]. Five spectra were randomly selected as endmembers from the U.S. Geological Survey Digital Spectral Library [43], with 188 bands each, as shown in Fig. 5(a), the 3000 × 2400 subpixel label field was designed according to the land-use and land-cover features, corresponding to the number of selected endmembers, as shown in Fig. 5(b). The original high-resolution hyperspectral imagery was reconstructed by combining the selected endmembers with the simulated subpixel label field. Then, the mean-value processing of the original high-resolution imagery yields the coarse-resolution hyperspectral imagery with the downsampling scales of 4, 5, and 6, respectively. So, the reconstruction scale factor was set to 4, 5, and 6, respectively. Then, the simulated coarse-resolution hyperspectral imagery are further degraded using the zero-mean i.i.d. Gaussian noise. In this experiment, different bands in the hyperspectral images were contaminated by noise with the same variance using the following formula.
where x ie represents the value of the e band of the ith hyperspectral pixel without noise. N is the number of pixels in the simulated hyperspectral imagery. SN R is the signal-to-noise ratio, the larger the SN R, the less noise the imagery. In this experiment, the SNR was set to 10. Finally, the simulated  hyperspectral images were solved by nonnegative least squares, respectively, and the abundance maps of sizes 750 × 600 × 5, 600 × 480 × 5, and 500 × 400 × 5 were obtained as the input data. Fig. 6 shows a visual comparison of the experimental results. The visual effect obtained by the SACRF method is smoother than that by other SPM methods because the former considers the spatial prior model in the energy function and is constrained by the abundance in the implicit form. Given the strict abundance constraints, many noise artifacts exist in the results of PSSPM, SASPM, RBFSPM, and SSSPM, which did not provide satisfactory visual effects. The local magnification within the red and black rectangles is a clear display of various kinds of SPM results. Compared with those of other SPM methods, the local magnification map of the proposed method and the true label [see Fig. 5(a)] are the most consistent. The noise artifact in the results of other methods seriously affects the visual effect. Table I shows the quantitative evaluation of the Kappa and OA of the SPM results for the true label. The OA and Kappa of the proposed SACRF method are higher than those of the other methods. The traditional PSSPM, SASPM, RBFSPM, and SSSPM methods are close because of their strict abundance constraints. The OA of the proposed method at the three scale factors are 97.27%, 97.08%, and 96.53%, which are 41.01%, 33.19%, and 33.05% higher than those of the SSSPM. Moreover, the OAs of the proposed method are 40.75%, 31.60%, and 32.47% higher than those of RBFSPM, 39.74%, 30.96%, and 31.17% higher than those of SASPM, and 41.61%, 33.44%, and 34.11% higher than those of PSSPM. Table II shows the producer accuracy statistics of experimental results from simulated data for various methods at different scale factors. The indicators reflected in the table are consistent with the results reflected in Table I and Fig. 6. It can be seen that the proposed method achieves good results and obtains high values for all classes of SPM. Table III shows the statistics of the running times of all algorithms to illustrate the operational efficiency of the algorithms for the simulated hyperspectral imagery. All algorithms in this study were run on matlab2010a, where the proposed method invoked the open-source toolbox GCO-3.0 written with C + + to implement the Graph-cut algorithm. As can be seen from Table III, the proposed method maintains operational efficiency despite not having the least running time when processing 3000 × 2400 sizes images under different scale factors.  Meanwhile, the proposed method can also handle the noise artifacts while maintaining high efficiency, which fully explains the advancement of the proposed method.

B. Experiment 2: Synthetic Hyperspectral Imagery
As shown in Fig. 7(a), the image of experiment 2 was taken in a small mountain village in Gulin, Sichuan Province, using a Cubert UHD185 hyperspectral camera. It includes three land cover classes: plastic film, soil, and vegetation. The hyperspectral camera is a frame camera, and the image consists of one high spatial resolution band of 1000 × 1000 and 138 low spatial resolution bands of 50 × 50. A hyperspectral image of 1000 × 1000 with 138 bands can be obtained through the internal fusion algorithm of the camera. The fused hyperspectral image was classified using the support vector machine algorithm in ENVI software. The wrong classification pixels were manually corrected to obtain an accurate 1000 × 1000 ground truth map and serve as the true label for the SPM, as shown in Fig. 7(c).
The region of interest was selected in the 1000 × 1000 hyperspectral image, as shown in Fig. 7(b). It was solved by nonnegative least squares to obtain 1000 × 1000 × 3 abundance maps. Then, the 1000 × 1000 × 3 abundance maps were downsampled at scale factors 2, 4, and 5. The 500 × 500 × 3, 250 × 250 × 3, and 200 × 200 × 3 synthetic abundance maps were obtained as input data. Therefore, the scale factors of the SPM algorithm in this experiment were set to 2, 4, and 5. Fig. 8 shows the experimental results for scale factor 4 of synthetic hyperspectral imagery. The results of PSSPM, SASPM, RBFSPM, and SSSPM methods strictly consider abundance constraints. Thus, many noise artifacts are produced, and the effect of SPM is seriously affected, such as the local magnification within the red and yellow rectangles, which is the result of the abundance error obtained by spectral unmixing. The SACRF method alleviates the effect of abundance error on the SPM by adding smoothed prior information, yielding smooth results, and showing desirable visual effects at different scale factors. Table IV lists the quantitative evaluation of these SPM methods at different scale factors to evaluate the effectiveness of these SPM methods using the OA and Kappa coefficients. Table IV shows that the traditional PSSPM, SASPM, RBFSPM, and SSSPM methods have similar quantitative evaluation indicators. The OA and Kappa of the proposed SACRF method are higher than those of the other methods. The OA of the proposed method  at the three scale factors are 88.45%, 88.62%, and 88.75%, which are 7.66%, 7.81%, and 7.73% higher than those of the SSSPM. Moreover, the OAs of the proposed method are 7.68%, 8.39%, and 8.54% higher than those of RBFSPM, 7.87%, 8.54%, and 8.78% higher than those of SASPM, and 8.38%, 10.23%, and 11.38% higher than those of PSSPM. Table V shows the producer accuracy of the experimental results of synthetic data for various methods at different scale factors. The indicators reflected in the table are consistent with the results reflected in Table IV and Fig. 8. The proposed method achieves good results and achieves high values for all classes. It can be noted that the producer accuracy of the vegetation is very low, because of its few pixels, but the proposed method still achieves the highest value.
The running times of all algorithms are counted in Table VI to illustrate the operational efficiency of the algorithms for the 1000 × 1000 imagery. Similar to experiment 1, experiment 2 was also run on matlab2010a, where the proposed method invoked the open-source toolbox GCO-3.0 written with C + + to implement the Graph-cut algorithm. Table VI shows that the proposed method has the shortest running time and the highest efficiency when processing the 1000 × 1000 sizes image at different scale factors. The effectiveness of the proposed method is fully illustrated.

C. Experiment 3: Real Multispectral Imagery
The data for this experiment were obtained using the multispectral images taken from the Landsat 8 satellite, which was developed by the United States and then launched on February 11, 2013. Landsat 8 satellite carries two main loads: Operational Land Imager (OLI) and Thermal Infrared Sensor. OLI has nine bands. The spatial resolution of band1−band7 and band9 is 30 m. Band8 is a panchromatic band, and its spatial resolution is 15 m. In this experiment, the data of the first seven bands of Landsat8 images were processed, and two regions were selected for experiments.
The data areas were divided into three classes: town, water, and vegetation. The multispectral band with a spatial resolution of 30 m and the panchromatic band with a spatial resolution of 15 m were fused using the Gram-Schmidt Pan Sharpening fusion   algorithm brought by ENVI software to obtain the true label and facilitate the quantitative evaluation of the SPM algorithm by obtaining the 920 × 1160 multispectral image with a spatial resolution of 15 m, as shown in Fig. 9(b). Then, the fused multispectral image was classified using the SVM algorithm in ENVI software. The wrong classification pixels were manually adjusted and finally obtained a true label, with a size of 920 × 1160, as shown in Fig. 9(c). Thus, the scale factor of this experiment was also set to 2. As shown in Fig. 10, this experiment selected the region of interest on the original multispectral image of 460 × 580 to obtain the input abundance. Then, unmixing was performed  using a nonnegative least square algorithm to obtain an abundance map measuring 460 × 580 × 3 as the input data, as shown in Fig. 10(b)-(d).
The results of the SPM in Fig. 11 are consistent with the results of experiments 1 and 2. The SPM result obtained by the SACRF method proposed in this experiment is smoother than that by other SPM methods. Given that PSSPM, SASPM, RBFSPM, and SSSPM methods strictly follow the abundance constraints, many noise artifacts exist in the results, seriously affecting the final SPM results. The local magnification within the red and black rectangles in Fig. 11 shows that the results obtained by the proposed SPM method best fit the true label. Moreover, the proposed SPM method eliminates many noise artifacts and has the best visual effect. Table VII shows the quantitative evaluation of the Kappa and OA. The proposed SPM method has the highest accuracy value and the best quantitative evaluation index. In particular, its OA is 86.33%, which is 7.54%, 11.39%, 11.05%, and 12.16% higher than the OA of SSSPM, RBFSPM, SASPM, and PSSPM, respectively.
Table VIII is a producer accuracy statistical table for the experimental results of region Baodi of each method at the scale factor of 2. The indicator reflected in the table is consistent with the results shown above in Table VII and Fig. 11. It can be seen that the proposed method achieves high producer accuracy for the SPM of all classes. Table IX shows the time used to process the image. Similar to the images in the above two experiments, the image in   part was also run on matlab2010a, where the proposed method invoked the open-source toolbox GCO-3.0 written with C + + to implement the Graph-cut algorithm. Table IX shows that the proposed SPM method has the shortest running time and is the most efficient for the image.
2) The region of Heze, Shandong, China: This region was acquired on September 2, 2016, in Heze, Shandong Province, and the size is 510 × 550, as shown in Fig. 12(a). The data were divided into three land cover classes: town, water, and vegetation. The multispectral band with a spatial resolution of 30 m and the panchromatic band with a spatial resolution of 15 m were fused using the Gram-Schmidt Pan Sharpening fusion algorithm brought by the ENVI software to obtain the true label and facilitate the quantitative evaluation of the SPM algorithm by obtaining a 1020 × 1100 size multispectral image with a spatial resolution of 15 m, as shown in Fig. 12(b). Then, the fused multispectral image was classified using the SVM algorithm in ENVI software. The wrong classification pixels was manually adjusted, finally obtaining a true label map measuring 1020 × 1100, as shown in Fig. 12(c). Thus, the scale factor for this part of the experiment is set to 2.
As shown in Fig. 13(a), this experiment selected the region of interest on the original multispectral image of 510 × 550 to obtain the input abundance. Then, unmixing was performed using a nonnegative least square to obtain an abundance map measuring 510 × 550 × 3 as the input data, as shown in Fig. 13(b)-(d). Fig. 14 shows a visual comparison of the experimental results. The visual effect obtained by the SACRF method is smoother than that by other SPM methods because the former considers the spatial prior model in the energy function and is constrained by the abundance in the implicit form. Given the strict abundance constraints, many noise artifacts exist in the results of PSSPM, SASPM, RBFSPM, and SSSPM, which did not provide satisfactory visual effects. The local magnification within the red and yellow rectangles is a clear display of various kinds of SPM results. Compared with those of other SPM methods, the local magnification map of the proposed method and the true label are the most consistent. The noise artifact in the results of other methods seriously affects the visual effect.     achieves high producer accuracy for the SPM of all classes, which illustrates the effectiveness of the proposed method.
The runtimes of all algorithms are counted in Table XII to illustrate the operational efficiency of the algorithms for the 1020 × 1100 image. As with the above experiments, this experiment was also run on matlab2010a, where the proposed method invoked the open-source toolbox GCO-3.0 written with C + + to implement the Graph-cut algorithm. Table XII shows that the proposed method has the shortest running time and the highest efficiency. The effectiveness of the proposed method is fully illustrated.

D. Sensitivity Analysis of λ
The tuning parameter λ in the proposed method is vital in controlling the weight between unary and pairwise potential terms. The smaller the λ value, the greater the spatial smoothing effect played by the pairwise term, and the smoother the SPM result, the larger the λ value, the greater the unary term plays, and the weight of all land cover class probabilities to which the subpixel belongs slowly increases, leading to a slow increase in noise artifacts. A hyperparameter search with interval 1 was performed in the range of 1-20 in all the datasets described above to determine the optimal λ values, and its OA was calculated. Fig. 15 shows the statistical summary results of the hyperparameter search.  As can be seen from Fig. 15, with the increasing λ value, the accuracy first increases and then decreases, and finally stabilizes, which meets ours expectation that the best balance point of the λ values between unary and pairwise potential terms can be found between 1-20.

E. Sensitivity Analysis of Noise
In order to examine the sensitivity of the different methods to the noise levels, we tested them on the simulated images with different noise levels measured with the SN R, i.e., SN R = 5, 10, 15, 20, and 25. The smaller the SN R value, the more noise in the imagery, and thus the greater the error of the abundance map produced by unmixing. In the experiments, we tested the results of SPM of different noise levels images at different scale factors, and the OA is shown in Fig. 16.
As can be seen from Fig. 16, the results of PSSPM, SASPM, RBFSPM, and SSSPM are not ideal when the image contains much noise. This is because the image contains more noise, resulting in a large error of the abundance map obtained from the unmixing, and thus the results of the SPM have low accuracy. This also shows that the dependence of these methods on the accuracy of the unmixing algorithms is great. The proposed method achieves good results at different noise levels, which shows that the proposed method can better handle the influence of the abundance error caused by unmixing on the SPM results.

V. DISCUSSION
Spectral unmixing technique estimates the area proportions of each land cover class within the mixed pixel (i.e., abundance).
SPM is a further technique for spectral unmixing, and under abundance constraints, the spatial distribution of each land cover class within the mixed pixel is obtained by considering spatial prior information, providing a land cover map with higher spatial resolution. Therefore, to satisfy the abundance constraint is an important problem in SPM, which is the key to maintain the significance of the SPM algorithm. Fig. 1 shows that there are usually two ways for SPM techniques to preserve the abundance constraint. One is a two-step strategy method, which requires spectral unmixing first to obtain the abundance information, and then performs the SPM processing under the abundance constraints of display form, that is, during the class allocation, the number of subpixels of each class in the mixed pixel is enforced according to the abundance information. The two-step strategy method has the advantage of high operational efficiency, but it is very vulnerable to abundance errors. The other is the optimization strategy method, which constructs the likelihood term in the objective function through a linear mixed model, thus transforming the abundance constraint in the SPM into a way of linear spectral reconstruction method, namely, minimize the difference between the product between abundance and endmembers and the observed mixed pixel spectrum. Although this method relate restriction on abundance and eliminates the uncertainty caused by the spectral unmixing algorithm, the operation efficiency is slow.
In order to make the SPM algorithm have the advantage of high operation efficiency, and can also eliminate the influence of abundance error, this article uses CRFs to model the SPM. However, when constructing the unary potential term in the energy function of the CRFs, the key lies in how to satisfy the abundance constraint while keeping the probability that subpixel belongs to each class. To meet this requirement, this article designs the SAAM. The traditional spatial attraction model quantifies the spatial correlation as a subpixel within the central mixed pixel subject to the spatial attraction of different land cover classes in the neighborhood mixed pixel. But the model belongs to the two-step strategy method, and the abundance constraint is in the display form [i.e., (5)], which enforces the number of subpixels of each class in the mixed pixel, and cannot be applied to the CRFs. The SAAM obtains the spatial adaptive attraction value by adaptive adjusting the spatial attraction value. The spatial adaptive attraction value adjusts the display form of the abundance constraints of spatial attraction value to an implicit form [i.e., (7)] with an adaptive scheme, while maintaining the relative size of the spatial correlation of each land cover class. Thus, it can be applied to the CRFs to meet the requirements of the SPM technology in constructing the unary potential term in the energy function of the CRFs.
In order to verify the efficiency advantages of the proposed method, the spatial size of the experimental data is relatively large, so the algorithm operation efficiency is verified by calculating the running time. However, when running on largesize data, the operation time of the optimization strategy SPM method is too long to make statistics, so the comparison method in this article belongs to the two-step strategy SPM method. Experiments show that the proposed algorithm guarantees the operational efficiency. Meanwhile, to illustrate the advantages of the proposed algorithm in eliminating the uncertainty caused by the spectral unmixing algorithm, this article takes the abundance of the nonnegative least squares solution as the input data, and the accuracy is not the highest, and the abundance is solved with obvious error. It can be concluded from the experiment that the comparison methods are seriously affected by the abundance error, so there are a lot of noise artifacts in the results, and the accuracy is low. And the proposed algorithm in this article can effectively solve the uncertainty caused by the unmixing algorithm when comparing it with other methods, thus achieving the results with relatively high accuracy than the comparison method.
On the other hand, our approach also has some limitations. In this article, although the SAAM is designed to meet the requirements of the SPM technology in constructing the unary potential term in the energy function of the CRFs. But the SAAM is based on the improvement of the traditional attraction mode. The traditional attraction model is to obtain the spatial correlation of subpixels belonging to each class through abundance, which leads to that if the abundance error is too large, it will affect the accuracy of the spatial correlation of subpixels belonging to each class, thus affecting the proposed method. Therefore, in future studies, we can design more advanced unary potential term in the energy function of the CRFs while satisfying the abundance constraints by considering more abundant spatial prior information or linear spectral mixing information. At the same time, the proposed method in this article has the tuning parameter λ, which controls the balance between the unary and pairwise potential terms. If too small, the result will be too smooth, and if too large, the noise cannot be effectively processed. At present, an optimal parameter λ is still obtained through parameter search, so in the future research, an algorithm can be designed to find the optimal parameter λ in an automatic way.
It is worth noting that the research and application of deep learning in the field of remote sensing is more and more mature and extensive. Recently, some paper [44] has used the supervised learning method to train the designed deep learning network by using labeled airborne data and unlabeled spaceborne data, so as to realize the mapping of urban land cover with high spatial resolution by using multisource data with different spatial resolution and spectral characteristics, providing a new perspective for large-scale land cover mapping. In the next research, we will further combine the theory of subpixel mapping with deep learning, and make breakthroughs in reducing training samples, improving generalization capability and improving accuracy.

VI. CONCLUSION
In this study, the SPM method based on the SAAM and CRFs are proposed. First, an SAAM is proposed. Spatial adaptive attraction values are obtained by adaptively adjusting the spatial attraction values obtained based on the conventional spatial attraction model. Spatial adaptive attraction values turn the display form of the abundance constraints into an implicit form for expression, allowing it to maintain the relative size of the spatial correlation of land cover class and characterize the limit of the abundance constraints. Thus, subsequent smoothing of noise artifacts is facilitated. Second, the spatial adaptive attraction values and spatially smoothed prior are modeled in CRFs. The unitary potential designed from the spatial adaptive attraction values can calculate the cost that each subpixel is given the corresponding class label. The binary potential modeled by a MLL model favors the adjacent subpixels to adopt the same class labels to achieve a smooth effect. Finally, the CRFs model is optimized using the Graph-cut, so that the proposed SPM can not only ensure the advantages of fast running speed and high efficiency of the two-step strategy SPM method, but also extinguish the influence of abundance error to a certain extent and decrease the noise artifact on the results of SPM. The experiments of the three images indicate that the result generated by the proposed SPM method is more reliable than the results generated by the four previous SRM methods. In future work, the design of unitary and binary potential in CRFs should be further explored for applications to remote sensing image processing in specific fields.