A New Spatial Attraction Model for Improving Subpixel Land Cover Classification

Subpixel mapping (SPM) is a technique that produces hard classification maps at a spatial resolution finer than that of the input images produced when handling mixed pixels. Existing spatial attraction model (SAM) techniques have been proven to be an effective SPM method. The techniques mostly differ in the way in which they compute the spatial attraction, for example, from the surrounding pixels in the subpixel/pixel spatial attraction model (SPSAM), from the subpixels within the surrounding pixels in the modified SPSAM (MSPSAM), or from the subpixels within the surrounding pixels and the touching subpixels within the central pixel in the mixed spatial attraction model (MSAM). However, they have a number of common defects, such as a lack of consideration of the attraction from subpixels within the central pixel and the unequal treatment of attraction from surrounding subpixels of the same distance. In order to overcome these defects, this study proposed an improved SAM (ISAM) for SPM. ISAM estimates the attraction value of the current subpixel at the center of a moving window from all subpixels within the window, and moves the window one subpixel per step. Experimental results from both Landsat and MODIS imagery have proven that ISAM, when compared with other SAMs, can improve SPM accuracies and is a more efficient SPM technique than MSPSAM and MSAM.


Introduction
Land use and land cover (LULC) information is very important in many scientific studies and applications.Remote sensing is the only feasible way of obtaining LULC information for large geographic areas.Many algorithms have been developed to classify various remote sensing data to obtain LULC maps [1][2][3][4].However, remote sensing images often contain mixed pixels, since the sensor's instantaneous field of view (IFOV) includes more than one land cover class [5,6].The existence of mixed pixels leads to three main problems which need be solved: (1) What classes of land cover does a mixed pixel contain?(2) What are the proportions of land cover classes in a pixel?(3) What is the subpixel spatial distribution of land cover classes [7]?For the first problem, end-member extraction algorithms [8], such as the pixel purity index (PPI), N-FINDER, iterative error analysis (IEA), etc., have been developed.For the second problem, soft or fuzzy classification algorithms, which allocate all classes, in varying proportions, to each pixel [9], have been proposed.These algorithms can be broadly categorized as linear spectral mixed models (LSMMs) and nonlinear mixture models [10][11][12].Due to the intrinsic complexity of the mixture modeling and the difficulty in obtaining scene parameters, nonlinear mixture models (NLMMs) have not been applied as widely as LSMMs [10,11].For the third problem, subpixel mapping (SPM) or super resolution mapping (SRM) methods, which produce hard classified maps at a spatial resolution finer than that of the input images [13], have been developed in recent decades.This paper focuses on SPM, which has been proved as an alternative method for obtaining land use/land cover with an acceptable accuracy [14][15][16][17][18][19].
In 1993, Schneider first introduced a knowledge-based analysis technique for the automatic localization of field boundaries in agricultural areas [20].Atkinson formally proposed the concept of SPM and mentioned that SPM can be considered as the post-processing of soft classification based on spatial dependence theory [21].With the assumption of spatial dependence, Verhoeye et al. [6] proposed a spatial dependence mathematical model and employed linear optimization techniques to find the maximum of the dependence.Considering each pixel as a neuron, Tatem et al. [22][23][24] and Wang et al. [25] applied a Hopfield neural network (HNN) to map subpixel land cover.The HNN increases the spatial correlation between neighboring subpixels and minimizes iterations to obtain SPM results.Atkinson [13] utilized a two-point histogram method to optimize the match between the target and the current realization of the two-point histograms for subpixel classes within pixels.Atkinson [9,26] developed a pixel swapping algorithm (PSA) to maximize the spatial correlation between neighboring subpixels, by changing the spatial arrangement of subpixels.PSA was initially designed to work for binary-class images, and was later expanded to work on multiple-class images [27][28][29].Mertens et al. [30] proposed a subpixel/ pixel spatial attraction model (SPSAM) to calculate the spatial attractions between subpixels and their neighboring pixels.In addition, the Markov random field [31,32], genetic algorithms [33,34], and indicator cokriging-based geostatistical methods [35,36], have also been successfully applied in SPM.
Among all of the aforementioned SPM methods, SPSAM has several advantages in terms of both its simplicity and its explicit physical meanings [30].For example, spatial attraction, which calculates the spatial correlation between subpixels and their surrounding pixels, is used in SPSAM as a simple tool to directly convey spatial dependence.Without requiring prior knowledge on the spatial structure, which is essential in some learning-based SPM methods, such as a two-point histogram, indicator cokriging-based methods, and genetic algorithms, etc., SPSAM can obtain satisfactory SPM results.However, SPSAM ignores the uncertainty of the spatial distribution of subpixels within surrounding pixels and fails to adequately consider the spatial correlation between subpixels within the central pixel.Consequently, the SPM results obtained by SPSAM are noisy and its accuracy is limited [37].Therefore, based on SPSAM, Wang et al. proposed a modified SPSAM (MSPSAM) [37] which estimates the spatial attractions according to the distribution of subpixels within neighboring pixels.They also proposed a mixed spatial attraction model (MSAM) for improving the SPM result.MSAM integrates the spatial attraction from both the immediate surrounding subpixels of the current subpixel within the central pixel, and all of the subpixels within the immediate neighboring pixels of the central pixel.
Among the abovementioned spatial attraction models (SAMs), MSAM is the only SAM which considers the spatial attractions generated by not only the subpixels of neighboring pixels, but also the subpixels of the central pixel.However, MSAM assumes that the spatial attraction of a subpixel being considered (called the current subpixel hereafter) is only influenced by the neighboring subpixels, instead of all the subpixels within the central pixel.Yet, according to the spatial dependence theory, the subpixels within the central pixel can exert stronger spatial attractions than those within neighboring pixels because they are spatially closer to the current subpixel.In MSAM, after spatial attraction of all subpixels within the central pixel is calculated, the calculation moves to the next pixels.As a result of this, the same set of subpixels within the neighboring pixels is used in calculating the spatial attraction of all the subpixels within the central pixel.When the current pixel, of which the spatial attraction is being calculated, is not located at the center of the central pixels, the subpixels within the pixel next to the neighboring pixels, which are closer to the current pixel than many of the subpixels within the neighboring pixels, are not considered.This results in the unequal treatment of attraction from surrounding subpixels of the same distance.
This study proposes an improved spatial attraction model (ISAM) to overcome the shortcomings of MSAM.ISAM estimates the attraction value of the current subpixel at the center of a moving window from all the subpixels within the window, and moves the window one subpixel per step, to improve the SPM result.

Methodology
2.1.Subpixel Mapping (SPM): Theory SPM aims to determine the most likely locations of the class fractions within a pixel.The general methods of SPM include three steps: (1) Utilizing spectral mixture analysis models to obtain soft class fraction (proportion) images at an original (coarse resolution) pixel resolution; (2) Dividing the original pixels into a series of subpixels, assuming that one subpixel only contains a specific class, to determine the number of subpixels for each class; (3) Applying spatial distribution features of classes and other prior knowledge, to map the subpixel spatial distribution of classes.From the abovementioned steps, it is obvious that the spatial distribution features of classes are the critical factor of SPM.A random subpixel distribution of classes can be assumed if prior knowledge is lacking.However, according to the spatial dependence theory, the land covers of two adjacent subpixels are more similar than those of two distant subpixels.Therefore, Atkinson considered the spatial dependence theory as the basis for SPM [21].
Figure 1 illustrates the spatial dependence theory of SPM.It shows a raster grid of 3 × 3 original (coarse resolution or mixed) pixels, with associated fractions of a specific class (Figure 1a).Each pixel is divided into S 2 subpixels (S is scale factor), each corresponding to 1/ S 2 area of the original pixel.Although both Figure 1b,c can present the possible results of the subpixel allocation of the gray class corresponding to the indicated proportion in Figure 1a, according to spatial dependence theory, this is more likely to represent the ground truth.This study proposes an improved spatial attraction model (ISAM) to overcome the shortcomings of MSAM.ISAM estimates the attraction value of the current subpixel at the center of a moving window from all the subpixels within the window, and moves the window one subpixel per step, to improve the SPM result.

Subpixel Mapping (SPM): Theory
SPM aims to determine the most likely locations of the class fractions within a pixel.The general methods of SPM include three steps: (1) Utilizing spectral mixture analysis models to obtain soft class fraction (proportion) images at an original (coarse resolution) pixel resolution; (2) Dividing the original pixels into a series of subpixels, assuming that one subpixel only contains a specific class, to determine the number of subpixels for each class; (3) Applying spatial distribution features of classes and other prior knowledge, to map the subpixel spatial distribution of classes.From the abovementioned steps, it is obvious that the spatial distribution features of classes are the critical factor of SPM.A random subpixel distribution of classes can be assumed if prior knowledge is lacking.However, according to the spatial dependence theory, the land covers of two adjacent subpixels are more similar than those of two distant subpixels.Therefore, Atkinson considered the spatial dependence theory as the basis for SPM [21].
Figure 1 illustrates the spatial dependence theory of SPM.It shows a raster grid of 3 × 3 original (coarse resolution or mixed) pixels, with associated fractions of a specific class (Figure 1a).Each pixel is divided into S 2 subpixels (S is scale factor), each corresponding to 1/ S 2 area of the original pixel.Although both Figure 1b,c can present the possible results of the subpixel allocation of the gray class corresponding to the indicated proportion in Figure 1a, according to spatial dependence theory, this is more likely to represent the ground truth.Based on Atkinson's study [21], Verhoeye et al. [6] presented a mathematical model for SPM, which transformed the SPM problem into one of assigning classes to the subpixels using linear optimization techniques.Suppose that the coarse resolution pixels are to be divided into S 2 subpixels.The number of subpixels that have to be assigned to class c is NSPc and has been derived from soft class fraction images.A measure for spatial dependence SDVc,j will be computed for class c at each subpixel j.Each subpixel has to be assigned a value of one or zero for each class, one indicating an assignment to a particular class.Following this, the problem of assigning each subpixel to a specific class emerges, which has the maximum value of the spatial dependence.
The mathematical model can be expressed as Equation (1): where ∈ 1,2, … , , s the total number of classes in the study case; , the model key parameter, respresents the spatial dependence values (SDV) of subpixel when it is assigned to Based on Atkinson's study [21], Verhoeye et al. [6] presented a mathematical model for SPM, which transformed the SPM problem into one of assigning classes to the subpixels using linear optimization techniques.Suppose that the coarse resolution pixels are to be divided into S 2 subpixels.The number of subpixels that have to be assigned to class c is NSP c and has been derived from soft class fraction images.A measure for spatial dependence SDV c,j will be computed for class c at each subpixel j.Each subpixel has to be assigned a value of one or zero for each class, one indicating an assignment to a particular class.Following this, the problem of assigning each subpixel to a specific class emerges, which has the maximum value of the spatial dependence.
The mathematical model can be expressed as Equation ( 1): where c ∈ {1, 2, . . . ,C}, C is the total number of classes in the study case; SDV cj , the model key parameter, respresents the spatial dependence values (SDV) of subpixel p j when it is assigned to class c; x cj , the choice variable of subpixel p j when it is assigned to class c, is defined in Equation ( 2): The model meets two constraints: Equation ( 3) means that only one subpixel can be assigned to a specific class and Equation ( 4) means that the number of subpixels belonging to class c within an original pixel has to be NSP c .

SPSAM
SPSAM was proposed based on the theory of spatial dependence.The assumptions of SPSAM are: (1) the fraction values of neighboring pixels exert the attraction toward subpixels within a central pixel; (2) a subpixel within the central pixel can only be attracted by pixels surrounding the central one; and (3) other pixels are assumed to be too distant to exert any attraction.Assuming that closer pixels attract the subpixel more than the distant ones, SPSAM calculates the attraction value of a neighboring pixel to a subpixel, based on their Euclidean distance.
Figure 2 illustrates the labeled pixels and subpixels, the coordinate system, and the distance between pixels and subpixels in SPSAM.A scale factor of S = 4 means that the original central pixel P 11 contains 16 subpixels with labels: p 44 , p 45 , p 46 , p 47 , . . ., p 74 , p 75 , p 76 , p 77 .The distance between subpixel p ij within the central pixel and the neighboring pixel P MN can be defined as: SDV c,ij , the attraction values of subpixel p ij , which is assigned class c, by the neighboring pixels, can be defined using Equation (6).
where F c (P MN ) is the fraction value of class c of the neighboring pixel P MN .
class c; , the choice variable of subpixel when it is assigned to class c, is defined in Equation ( 2 The model meets two constraints: Equation ( 3) means that only one subpixel can be assigned to a specific class and Equation ( 4) means that the number of subpixels belonging to class c within an original pixel has to be .

SPSAM
SPSAM was proposed based on the theory of spatial dependence.The assumptions of SPSAM are: (1) the fraction values of neighboring pixels exert the attraction toward subpixels within a central pixel; (2) a subpixel within the central pixel can only be attracted by pixels surrounding the central one; and (3) other pixels are assumed to be too distant to exert any attraction.Assuming that closer pixels attract the subpixel more than the distant ones, SPSAM calculates the attraction value of a neighboring pixel to a subpixel, based on their Euclidean distance.
Figure 2 illustrates the labeled pixels and subpixels, the coordinate system, and the distance between pixels and subpixels in SPSAM.A scale factor of S = 4 means that the original central pixel contains 16 subpixels with labels: , , , the attraction values of subpixel , which is assigned class c, by the neighboring pixels, can be defined using Equation (6).
where ( ) is the fraction value of class c of the neighboring pixel .

MSPSAM and MSAM
SPSAM estimates the spatial attractions of subpixel p ij by the neighboring pixel, ignoring that the other subpixels within the central pixel can exhibit spatial attraction toward the subpixel p ij , and the different spatial distribution of classes of the neighboring pixels can also have different spatial attraction values.Therefore, Wang et al., 2012b, proposed MSPSAM and MSAM to improve the SPM results.As illustrated in Figure 3a, MSPSAM calculates the spatial attractions of subpixel p ij within the central pixel for a given class by using all of the subpixels within the neighboring pixels that have been assigned to the same class as the given class of subpixel p ij .MSAM, illustrated in Figure 3b, computes the spatial attraction of subpixel p ij by using not only the subpixels within the neighboring pixels, but also the neighboring subpixels of subpixel p ij within the central pixel.3a, MSPSAM calculates the spatial attractions of subpixel within the central pixel for a given class by using all of the subpixels within the neighboring pixels that have been assigned to the same class as the given class of subpixel .MSAM, illustrated in Figure 3b, computes the spatial attraction of subpixel by using not only the subpixels within the neighboring pixels, but also the neighboring subpixels of subpixel within the central pixel.

Improved Spatial Attraction Model (ISAM)
Referring to PSA in Atkinson [9,26] and the modified pixel swapping algorithm (MPSA) in [29], ISAM computes the spatial attractions of the current subpixel at the center of a moving window by the surrounding subpixels within the window, which can be 2S + 1 times the size of the subpixel and moves the window one subpixel per step.As an illustration of ISAM, Figure 4 shows that subpixel within the central pixel is attracted by the surrounding subpixels within a moving window that is double the size of the original pixel.In Figure 4, S is the scale factor, which divides an original pixel into S 2 subpixels, each corresponding to 1/S 2 area of the original one.

Improved Spatial Attraction Model (ISAM)
Referring to PSA in Atkinson [9,26] and the modified pixel swapping algorithm (MPSA) in [29], ISAM computes the spatial attractions of the current subpixel at the center of a moving window by the surrounding subpixels within the window, which can be 2S + 1 times the size of the subpixel and moves the window one subpixel per step.As an illustration of ISAM, Figure 4 shows that subpixel p ij within the central pixel P cen is attracted by the surrounding subpixels within a moving window that is double the size of the original pixel.In Figure 4, S is the scale factor, which divides an original pixel into S 2 subpixels, each corresponding to 1/S 2 area of the original one.

Improved Spatial Attraction Model (ISAM)
Referring to PSA in Atkinson [9,26] and the modified pixel swapping algorithm (MPSA) in [29], ISAM computes the spatial attractions of the current subpixel at the center of a moving window by the surrounding subpixels within the window, which can be 2S + 1 times the size of the subpixel and moves the window one subpixel per step.As an illustration of ISAM, Figure 4 shows that subpixel within the central pixel is attracted by the surrounding subpixels within a moving window that is double the size of the original pixel.In Figure 4, S is the scale factor, which divides an original pixel into S 2 subpixels, each corresponding to 1/S 2 area of the original one.Therefore, defined in Equation ( 7), J c,ij , the spatial attraction of subpixel p ij which is assigned to class c (c can be any class in the study case), can be estimated as an inverse-distance weighted function of its surrounding subpixels within the moving window.
where d p mn , p ij is the distance between subpixel p ij and the surrounding subpixel p mn defined in Equation ( 8), and x c,ij is a choice variable defined in Equation (9).
x c,ij = 1, i f subpixel p ij and p mn are assigned to the same class c 0, otherwise The mathematical model thus becomes Equation ( 10): The model constraints of ISAM are shown as Equations ( 11) and ( 12): Equation ( 11) means that only one subpixel can be assigned a specific class, and Equation ( 12) means that the number of subpixels that have to be assigned class c is n c , which has been determined from the fraction images.
Based on spatial dependence theory, maximum J will be retrieved when the top n c subpixels within the central pixel are assigned to class c.
The ISAM proposed in this study is similar to the MPSA proposed by Shen et al., 2009 [29], but there are three differences between them, including: (1) the initialization of subpixel class allocation.The former uses the random allocation and the latter allocates the subpixel class based on spatial attractiveness; (2) the size of the moving window.The former's is 2S + 1 times the size of the sub-pixel, and the latter's is unfixed; and (3) the calculation of the distance weighting parameter.The former takes the Euclidean distance as the weighting parameter of two subpixels and the latter uses an exponential model containing the Euclidean distance and a non-linear parameter.These differences make ISAM computationally more efficient.
As seen for other SPM methods, theoretically, ISAM also has a few limitations, including: (1) the classification accuracy of ISAM depends on the accuracy of the soft class fraction image, which is the result of spectral mixture analysis; and (2) ISAM prefers a subpixel mapping H-resolution-case image, of which the pixels are smaller than the objects of interest or landcover, to the L-resolution-case one, of which the pixels are much larger than the objects of interest or landcover [38,39].

The Algorithm of ISAM
The algorithm flowchart of ISAM is shown in Figure 5, and the detailed steps of the algorithm are given below:

Experiments and Results
To test and validate the advantages of the proposed ISAM algorithm, this study compared the SPM accuracies of ISAM with those of other SAMs by using Landsat OLI and MODIS imagery.ISAM, SPSAM, MSPSAM, and MSAM were implemented with ENVI IDL 8.3, and the following experiments were all accomplished in ENVI 5.1.A workstation computer with an Intel Quad 2.67 GHz processor and 4 GB RAM was used for this study.In all experiments, the moving window for ISAM is set to two original pixels plus one subpixel in size, so that the current subpixel is always at the center of the window.

Data Sets
A scene of the Landsat-8 Operational Land Imager (OLI) image with its identifier LC81460322015267LGN00 was downloaded from the U.S. Geological Survey (USGS) official website [40] and the region of interest (ROI) of 552 × 424 pixels in size was subset from the scene.By applying a support vector machine (SVM, of which the training samples were selected through visually interpreting both the high spatial resolution Google earth and the Landsat-8 images) to the Landsat image of ROI, a hard land cover classification map, which is used as a reference map, was created for the ROI at 30-m resolution.Then, by using a bi-cubic resampling algorithm, the land cover map at the original 30-m resolution was aggregated to 60-, 120-and 240-m resolution to form soft land cover fraction maps which were 276 × 212, 138 × 106, and 69 × 53 pixels in size, respectively.

Experiment Results
The accuracy is measured by comparing the SPM results from ISAM, SPSAM, MSPSAM, and MSAM with the reference map.The Overall Accuracy (OA) and Kappa Coefficient (κ), the most commonly used indices for classification accuracy assessment, are applied in this study.The experimental results are displayed in Figures 6-8 and the accuracy measures are shown in Tables 1-3.From Figures 6-8, a visual assessment of the image quality reveals that: (1) the SPM of the four SAM algorithms can reveal more details than hard classification results at scale factors of S = 2, 4, and 8; (2) the SPM from ISAM can obtain more details than those from other SAMs (see the details within red rectangles in Figure 8c).From Tables 1-3, a quantitative comparison analysis proves that: (1) the OA values of ISAM and the other SAMs at different scale factors are all above 95%, while all κ values are greater than 0.92, which means that all SAMs are effective SPM techniques; (2) both the OA and κ values of ISAM are the highest among all SAMs; (3) both the OA and κ values of all SAMs decrease with an increase of the scale factor, since the bigger the scale factor is, the coarser the aggregated image can be, and the less the detail of an initial image can convey; (4) the OA and κ of SPSAM are lower than those of MSPSAM and MSAM when the scale factor, S, is 2 or 4, while OA and κ of SPSAM are higher than those of MSPSAM and MSAM when the scale factor, S, is 8. From both visual and quantitative assessments, it is concluded that ISAM is a more accurate SPM technique than the other SAMs.From Tables 1-3, it can be also found that OA and κ of the four SAMs decreases by about 1.816%~4.035%and 0.029~0.065,respectively, when the scale factor S increases from 2 to 8. The sensitivity of the scale factor to the accuracy of the four SAMs is not obvious in this experiment and S can be determined according to the goal of experiment.The comparison results of computational efficiency among ISAM, SPSAM, MSPSAM, and MSAM are shown in Table 4 when the scale factor S is set to 8 and the image size is 424 × 552 subpixels.Among the four SAMs, SPSAM needs the fewest number of iterations and the least time to achieve the optimal results, while ISAM needs three more iterations than SPAM, but one and five iterations fewer than MSPSAM and MSPSAM, respectively.Furthermore, ISAM reduces to almost 900 s per iteration or 2.2 h and 7.8 h less per entire computation than MSPSAM and MSPSAM.In other words, ISAM is more efficient than MSPSAM and MSPSAM, but less efficient than SPSAM.The reason for this is that SPSAM only considers the attraction from eight neighboring pixels at the original resolution to the subpixels of the central pixels, resulting in a significant reduction in the computational requirement, while the others all consider the attraction from the subpixels of the neighboring pixels.Although the computing power is not a significant limit factor nowadays for processing remote sensing images, timesaving is significant if applying the ISAM algorithm to large-scale subpixel hard classification, which may need to classify hundreds of remote sensing images, instead of MSPSAM and MSAM.Nevertheless, if the computational efficiency is the main concern, SPSAM is the better choice since it is 67 times faster than ISAM.The experiment in this section, taking the land cover classification results from the Landsat OLI image with the support vector machine (SVM) as the reference map, focuses on extracting the subpixel land cover from MODIS imagery.
The MODIS image, which covered the area of the Landsat OLI image and was acquired on 16 May 2014, was downloaded from the U.S. NASA official website [41].The MODIS image was re-projected from Sinusoidal to Universal Transverse Mercator (UTM) projection with the MODIS Re-projection Tool (MRT), so that it could be co-registered with the Landsat image.Smoothing filter-based intensity modulation (SFIM) [42] was utilized to sharpen the MOD09GA bands 3-7 data to a 250-m pixel resolution of the MOD09GQ bands 1-2.After using the pixel purity index (PPI) method to select the end member of the land cover classes, the linear spectral mixing model (LSMM) was applied to extract the soft classification fraction from MODIS bands 1-2 and sharpened bands 3-7 data.
The Landsat OLI image (data identifier: LC81460322014136LGN00), which was acquired on the same day as the MODIS image, was downloaded from the U.S. Geological Survey (USGS) official website [40].The Landsat image then was strip-repaired, atmospherically corrected, ROI clipped, and resampled into a pixel resolution of 31.25 m, to meet the requirement of SPM from the 250-m resolution MODIS image with an integral number of scale factors.The preprocessed image was used to extract the land cover classification map, which was used as the ground truth in this MODIS experiment, by utilizing the support vector machine (SVM).

Experiment Results
Setting scale factor S to 8, the subpixel resolution of the MODIS image will be 31.25 m, the same as the reference map.The soft classification fraction image is used as the input to the SPM algorithms of ISAM, SPSAM, MSPSAM, and MSAM.Similar to the experiment with the Landsat OLI image, the accuracy is measured by comparing the SPM results of ISAM and the other SAMs with the reference map.The SPMs from different SAM algorithms are displayed in Figure 9 and the accuracy measures are shown in Table 5.   9, it can be found that the SPM map from ISAM contains more details than those from other SAMs.For instance, the details within the red rectangular areas in the SPM Map from ISAM match the same areas in the reference map much better than those from other SAMs.Table 5 shows that the OA of the SPM results from ISAM is 82.44%, which is higher than those of other SAMs; and κ is 0.66, which is higher than that of SPSAM and is equal to those of MSPSAM and MSAM.From both a visual assessment and quantitative analysis, it can be concluded that ISAM is more effective than other SAMs in SPM.
The comparison results of computational efficiency among ISAM, SPSAM, MSPSAM, and MSAM are shown in Table 6.Among the four SAMs, SPSAM is the most efficient SPM technique, while ISAM is the second most efficient one.
Comparing ISAM with SPAM, OA and κ increase by 12.22% and 0.2, respectively, in the experiment with MOIDS (though the accuracies of ISAM increase slightly, compared with SPSAM, and MSPSAM and MSAM in the experiment with Landsat).Therefore, ISAM can improve SPM accuracies more than SPSAM and is a more efficient SPM technique than MSPSAM and MSAM.By making a visual assessment of SPM maps in Figure 9, it can be found that the SPM map from ISAM contains more details than those from other SAMs.For instance, the details within the red rectangular areas in the SPM Map from ISAM match the same areas in the reference map much better than those from other SAMs.Table 5 shows that the OA of the SPM results from ISAM is 82.44%, which is higher than those of other SAMs; and κ is 0.66, which is higher than that of SPSAM and is equal to those of MSPSAM and MSAM.From both a visual assessment and quantitative analysis, it can be concluded that ISAM is more effective than other SAMs in SPM.
The comparison results of computational efficiency among ISAM, SPSAM, MSPSAM, and MSAM are shown in Table 6.Among the four SAMs, SPSAM is the most efficient SPM technique, while ISAM is the second most efficient one.
Comparing ISAM with SPAM, OA and κ increase by 12.22% and 0.2, respectively, in the experiment with MOIDS (though the accuracies of ISAM increase slightly, compared with SPSAM, and MSPSAM and MSAM in the experiment with Landsat).Therefore, ISAM can improve SPM accuracies more than SPSAM and is a more efficient SPM technique than MSPSAM and MSAM.

Conclusions
In order to overcome the defects in the existing SAM techniques, this study has proposed an improved SAM (ISAM) for SPM, through extending the existing SAM techniques.Instead of computing the spatial attraction by the surrounding pixels in SPSAM, or by the subpixels within the surrounding pixels and the touching subpixels within the central pixel in MSAM, ISAM estimates the attraction of the current subpixel at the center of a moving window by using all of the subpixels within the window and moves the window one subpixel at a time.The design of the algorithm is more straightforward and logically consistent than existing SAM algorithms, resulting in the simplification of the algorithm implementation and efficiency in the algorithm execution.Experimental results from both Landsat and MODIS imagery show that ISAM improves the SMP accuracy over the existing SAMs and is computationally more efficient than SPSAM and MSAM.Overall, ISAM is an effective and efficient technique for SPM.

Figure 1 .
Figure 1.Illustration of spatial dependence theory for subpixel mapping in an 8 × 8 subpixel scene.(a) 3 × 3 coarse resolution pixels with the indicated proportion of a specific (gray) class; (b) and (c) The possible results of the subpixel allocation of the gray specific class.

Figure 1 .
Figure 1.Illustration of spatial dependence theory for subpixel mapping in an 8 × 8 subpixel scene.(a) 3 × 3 coarse resolution pixels with the indicated proportion of a specific (gray) class; (b) and (c) The possible results of the subpixel allocation of the gray specific class.

Figure 4 .
Figure 4. Illustration of the Improved Spatial Attraction Model (ISAM) (the size of the moving window is equal to 2S + 1).

Figure 4 .
Figure 4. Illustration of the Improved Spatial Attraction Model (ISAM) (the size of the moving window is equal to 2S + 1).

Figure 4 .
Figure 4. Illustration of the Improved Spatial Attraction Model (ISAM) (the size of the moving window is equal to 2S + 1).

Figure 5 .
Figure 5.The algorithm flowchart of the improved spatial attraction model (ISAM): Hmax-maximum steps of iteration; H-current step of iteration; Acc-the accuracies of SPM of previous iteration step; Acc_c-the accuracies of SPM of current iteration step; -the central pixel; -the subpixel which spatial attractions are currently calculated; C-the number of classes in the study case; Jc,ij-the spatial attractions of subpixel when it is assigned to class c, SLSc-the spatial location sequence of class c of pixel .

Figure 5 .
Figure 5.The algorithm flowchart of the improved spatial attraction model (ISAM): H max -maximum steps of iteration; H-current step of iteration; Acc-the accuracies of SPM of previous iteration step; Acc_c-the accuracies of SPM of current iteration step; P cen -the central pixel; p ij -the subpixel which spatial attractions are currently calculated; C-the number of classes in the study case; J c,ij -the spatial attractions of subpixel p ij when it is assigned to class c, SLS c -the spatial location sequence of class c of pixel P cen .

Figure 6 .
Figure 6.The comparison of SPM results among ISAM, SPSAM, MSPSAM, and MSAM (scale factor S = 2): (a) The classification result from Landsat data; (b) The hard classification result at scale factor S; (c-f) the results from ISAM, SPSAM, MSPSAM, and MSAM, respectively.

Figure 7 .
Figure 7.The comparison of SPM results among ISAM, SPSAM, MSPSAM, and MSAM (scale factor S = 4): (a) The classification result from Landsat data; (b) The hard classification result at scale factor S; (c-f) the results from ISAM, SPSAM, MSPSAM, and MSAM, respectively.

Figure 7 .
Figure 7.The comparison of SPM results among ISAM, SPSAM, MSPSAM, and MSAM (scale factor S = 4): (a) The classification result from Landsat data; (b) The hard classification result at scale factor S; (c-f) the results from ISAM, SPSAM, MSPSAM, and MSAM, respectively.

Figure 7 .
Figure 7.The comparison of SPM results among ISAM, SPSAM, MSPSAM, and MSAM (scale factor S = 4): (a) The classification result from Landsat data; (b) The hard classification result at scale factor S; (c-f) the results from ISAM, SPSAM, MSPSAM, and MSAM, respectively.

Figure 8 .
Figure 8.The comparison of SPM results among ISAM, SPSAM, MSPSAM, and MSAM (scale factor S = 8): (a) The classification result from Landsat data; (b) The hard classification result at scale factor S, (c-f) the results from ISAM, SPSAM, MSPSAM, and MSAM respectively.

Figure 8 .
Figure 8.The comparison of SPM results among ISAM, SPSAM, MSPSAM, and MSAM (scale factor S = 8): (a) The classification result from Landsat data; (b) The hard classification result at scale factor S, (c-f) the results from ISAM, SPSAM, MSPSAM, and MSAM respectively.

Figure 9 .
Figure 9.The comparison of subpixel mapping results among ISAM, SPSAM, MSPSAM, and MSAM (Scale factor S = 8; Mask represents the class of nonagricultural land; PML represents the class of plastic mulched landcover which is a type of farmland covered by plastic mulch film): (a) The classification results from Landsat data; (b) The hard classification results at scale factor S; (c-f) the results from ISAM, SPSAM, MSPSAM, and MSAM, respectively.

Figure 9 .
Figure 9.The comparison of subpixel results among ISAM, SPSAM, MSPSAM, and MSAM (Scale factor S = 8; Mask represents the class of nonagricultural land; PML represents the class of plastic mulched landcover which is a type of farmland covered by plastic mulch film): (a) The classification results from Landsat data; (b) The hard classification results at scale factor S; (c-f) the results from ISAM, SPSAM, MSPSAM, and MSAM, respectively.
SPSAM estimates the spatial attractions of subpixel by the neighboring pixel, ignoring that the other subpixels within the central pixel can exhibit spatial attraction toward the subpixel , and the different spatial distribution of classes of the neighboring pixels can also have different spatial attraction values.Therefore, Wang et al., 2012b, proposed MSPSAM and MSAM to improve the SPM results.As illustrated in Figure Sort SLSc (spatial location sequence for class c) by Jc,ij descending order Choose top subpixels and assign them to class c, according to SLSc Reassign the subpixels, which have been assigned to more than one class, to a specific class of which the spatial attractions reach the maximum.Make sure every subpixel is uniquely assigned to a specific class.
Set scale factor S Set maximum iteration step Hmax Randomly allocate subpixel to classes based on the pixel-level class fraction Estimate Acc, the accuracies of SPM of previous iteration step Initialize current iteration step H FOR each iteration//H FOR each pixel//Select a pixel as the central pixel Pcen