Combining Pixel Swapping and Simulated Annealing for Land Cover Mapping

Mixed pixels commonly exist in low-resolution remote sensing images, and they are the key factors hindering the classification of land covers and high-precision mapping. To obtain the spatial information at the subpixel level, subpixel mapping (SPM) technologies, including the pixel-swapping algorithm (PSA), that use the unmixed proportions of various land covers and allocate subpixel land covers have been proposed. However, the PSA often falls into a local optimum solution. In this paper, we propose a SPM method, the PSA_MSA algorithm, that combines the PSA and the modified simulated annealing algorithm to find the global optimum solution. The modified simulated annealing algorithm swaps subpixels within a certain range to escape the local optimum solution. The method also optimizes all the mixed pixels in a randomized sequence to further improve the mapping accuracy. The experimental results demonstrate that the proposed PSA_MSA algorithm outperforms the existing PSA-based algorithms for SPM. The images with different spatial dependences are tested and the results show that the proposed algorithm is more suitable for images with high spatial autocorrelation. In addition, the effect of proportion error is analyzed by adding it in the experiments. The result shows that a higher proportion error rate leads to larger degradation of the subpixel mapping accuracy. Finally, the performance of PSA_MSA algorithm with different ranges of selection on subpixel-swapping is analyzed.


Introduction
Remote sensing technologies have been increasingly used in many areas, such as environmental monitoring, agricultural production, and resource exploration. The spatial resolution and spectral resolution of spectral imaging systems are mutually restricted. To obtain a higher spectral resolution, the spatial resolution of a spectral imaging system has to be sacrificed [1]. As a result, the spatial resolution of spectral imaging systems is lower than that of panchromatic imagers. In images with low spatial resolutions, one image pixel might include several different types of land cover [2][3][4][5][6]. These land covers have different spectral characteristics. A pixel containing only one type of land cover is called a pure pixel. Otherwise, it is called a mixed pixel. To obtain the end-member information, spectral unmixing algorithms that can determine the proportion of different land covers in the mixed pixels have been proposed, including the linear mixture model [7,8], the fuzzy c-means classification model [9,10], and the neural network model [11,12]. Although proportional images can show the proportions of different land cover types in mixed pixels and provide a better understanding of the features contained in the mixed pixels, they still cannot provide the spatial distributions of the determined classes.

Principles of the PSA
We begin with a brief introduction of the principle behind and the functions of the PSA. The PSA constantly exchanges classes of subpixels to ensure that the subpixels from the same object have the largest spatial correlation. The PSA consists of three steps. Some symbols are shown in Figure 1.
Step 1. Calculate the gravitational value which reflects the spatial correlation of a mixed pixel. For a subpixel p i,j located at (i, j), its gravitational value can be calculated as follows: where A p i,j is the gravitational value of p i,j , N is the number of neighboring subpixels, λ k is the weighting function of distance, and z(p i,j , p k ) is the And operation function. The total gravitational value of a mixed pixel P a,b can be obtained by summing the gravitational values of all its subpixels p i,j .
Sensors 2020, 20, 1503 3 of 20 The total gravitational value of a mixed pixel Pa,b can be obtained by summing the gravitational values of all its subpixels pi,j. Step 2. Exchange the classes of the subpixels within the mixed pixel. For all subpixels of category 1, the subpixels with the lowest gravitational value are set to 0; for all subpixels of category 0, the subpixels with the lowest gravitational value are set to 1.
Step 3. Based on Equation (5), calculate the new total gravitational value of the mixed pixel. If the new gravitational value increases, then a subpixel-swapping result is retained. Otherwise, the locations of different land cover subpixels will return to the previous state.
Repeat step 1 to step 3 until the total gravitational value (Equation (5)) of the mixed pixel does not increase. The corresponding SPM results are considered to be the best localization results.
The current swapping method of the PSA is relatively simple, yet it is likely to fall into a local optimum solution of the objective function. To solve this problem, we introduce a SA algorithm to obtain the global optimum solution during subpixel swapping.

PSA with Simulated Annealing (SA)
The SA algorithm, which is used for solving optimization problems, was first proposed by Kirkpatrick et al. [35]. It is based on a strong analogy between the issue of solving optimization problems and the physical annealing of solids, which leads a system to its lowest energy state. Thus, the SA algorithm has the ability to escape from the local optimum solution [36,37].
A simulated annealing process includes the following steps: (1) select an objective function and set the starting temperature (simulated melting point), (2) simultaneously reduce the temperature and decrease the energy of the objective function, and (3) repeat step 2) until the target temperature reaches the ideal temperature and the objective function converges to the minimum value. Peng 32 directly applied the conventional SA algorithm to the swapping process of the PSA, and we refer to this algorithm as the PSA_SA. Peng used the spatial correlation to calculate the objective function used in the SA algorithm. The spatial correlation coefficients are calculated as follows: Step 2. Exchange the classes of the subpixels within the mixed pixel. For all subpixels of category 1, the subpixels with the lowest gravitational value are set to 0; for all subpixels of category 0, the subpixels with the lowest gravitational value are set to 1.
Step 3. Based on Equation (5), calculate the new total gravitational value of the mixed pixel. If the new gravitational value increases, then a subpixel-swapping result is retained. Otherwise, the locations of different land cover subpixels will return to the previous state.
Repeat step 1 to step 3 until the total gravitational value (Equation (5)) of the mixed pixel does not increase. The corresponding SPM results are considered to be the best localization results.
The current swapping method of the PSA is relatively simple, yet it is likely to fall into a local optimum solution of the objective function. To solve this problem, we introduce a SA algorithm to obtain the global optimum solution during subpixel swapping.

PSA with Simulated Annealing (SA)
The SA algorithm, which is used for solving optimization problems, was first proposed by Kirkpatrick et al. [35]. It is based on a strong analogy between the issue of solving optimization problems and the physical annealing of solids, which leads a system to its lowest energy state. Thus, the SA algorithm has the ability to escape from the local optimum solution [36,37].
A simulated annealing process includes the following steps: (1) select an objective function and set the starting temperature (simulated melting point), (2) simultaneously reduce the temperature and decrease the energy of the objective function, and (3) repeat step (2) until the target temperature reaches the ideal temperature and the objective function converges to the minimum value. Peng 32 directly applied the conventional SA algorithm to the swapping process of the PSA, and we refer to this algorithm as the PSA_SA. Peng used the spatial correlation to calculate the objective function used in the SA algorithm. The spatial correlation coefficients are calculated as follows: where N is the number of neighboring subpixels of each subpixel, while δ p i,j , p k is the And operation function. The final target optimization function of one mixed pixel is the summation of the spatial correlation coefficient of each subpixel within this mixed pixel and is defined as: The subpixel swapping process of PSA_SA is as follows: Randomly exchange the classes of subpixels in a mixed pixel, and calculate the new target optimization function A P a,b and the corresponding ∆A P a,b =A P a,b − A P a,b . If ∆A P a,b ≥ 0, the exchanged positions of subpixel land covers would be retained, and A P a,b = A P a,b . If ∆A P a,b < 0, calculate the acceptance probability ε = exp(−A P a,b /T) at the current SA temperature T, and compare ε with a random number r, which is in the range of (0,1). If ε > r, then the exchanged result would also be retained. Otherwise, the location of subpixels will return to the previous state.
Unlike the PSA, the PSA_SA contains an added SA algorithm, which has shown a great tolerance to local optima convergence and is often called a global optimizer that accepts deteriorated solutions with a certain probability in subpixel swapping [38]. Thus, the PSA_SA should be able to help mixed pixels to get rid of the problem of the local minimum in the process of increasing their gravitational value through subpixel swapping. However, the constant random exchange in the SA algorithm makes it difficult to use the subpixels of various land cover types to establish a reasonable spatial distribution. Therefore, the accuracy of the PSA_SA is not ideal, which will be demonstrated in the experimental results in Section 3.

PSA Based on the Modified Simulated Annealing Algorithm (PSA_MSA)
According to Atkinson, the purpose of the PSA is to maximize the spatial dependence, which is "the likelihood that observations close together are more alike than those that are further apart" 25, for subpixels of different land covers in a mixed pixel. To achieve this goal, his proposed PSA swaps the subpixels corresponding to the two types of land covers that have the smallest gravitational values. Peng's method merely swaps two randomly selected subpixels from all the subpixels of different land covers. Although this method considers the spatial dependence among the subpixels, the swapping range in this method is too wide and therefore cannot obtain sufficiently accurate results within a limited number of interactions.
In this research, based on the spatial dependence principle, we propose a PSA based on the modified simulated annealing algorithm. Figure 2a shows a mixed pixel with two kinds of land cover. Several subpixels with small gravitational values in the corresponding land cover form a group, and two groups of these subpixels are shown as slashed blocks in Figure 2b. In the proposed method, the SA algorithm randomly selects only one pixel from each group and swaps their classes. Figure 2c,d show two possible swapping results. This method constrains the subpixels, which are randomly swapped by SA, to subpixels with small gravitational values within their own classes. The procedure of the proposed PSA_MSA algorithm following various steps: Figure 3 illustrates the process of optimizing mixed pixels from step 1 to step 3 in the form of a flow chart. Step 1. Initialization According to the proportions of land cover in the mixed pixels, the subpixels with corresponding classes are randomly allocated in each mixed pixel. Calculate the total attractiveness of each mixed pixel by using either Equation (2) or Equation (4). Equation (2) considers the distance between the central subpixel and neighboring subpixels as the weight function, while Equation (4) treats The procedure of the proposed PSA_MSA algorithm following various steps: Figure 3 illustrates the process of optimizing mixed pixels from step 1 to step 3 in the form of a flow chart. The procedure of the proposed PSA_MSA algorithm following various steps: Figure 3 illustrates the process of optimizing mixed pixels from step 1 to step 3 in the form of a flow chart. Step 1. Initialization According to the proportions of land cover in the mixed pixels, the subpixels with corresponding classes are randomly allocated in each mixed pixel. Calculate the total attractiveness of each mixed pixel by using either Equation (2) or Equation (4). Equation (2) considers the distance between the central subpixel and neighboring subpixels as the weight function, while Equation (4) treats Step 1. Initialization According to the proportions of land cover in the mixed pixels, the subpixels with corresponding classes are randomly allocated in each mixed pixel. Calculate the total attractiveness of each mixed pixel by using either Equation (2) or Equation (4). Equation (2) considers the distance between the central subpixel and neighboring subpixels as the weight function, while Equation (4) treats neighboring subpixels with the same weighting value. We denote the algorithms based on Equations (2) and (4) as PSA_MSA1 and PSA_MSA2, respectively.
In the SA algorithm, set the values of the parameters such as the attenuation factor (a), the number of iterations (iter), the initial temperature T start , and the final temperature T stop . The cooling function is defined by T q+1 = T q × a, where q is the number of cooling times which reflects the number of iterations and T 0 = T start .
Step 2. Use the SA algorithm to randomly swap the selected subpixels According to the principle of spatial relevance, we can only swap subpixels with low attractiveness. Therefore, a constraint is applied to the SA algorithm by selecting subpixels that need to be swapped with attractiveness values within a certain range.
The low attractiveness range of is defined as follows: First, the subpixel gravitational values of two opposing classes are calculated separately. Here, we use M and N to represent two opposing classes. There are l subpixels in class M and S × S − l subpixels in class N.
Then, the gravitational values of M-class subpixels are sorted from smallest to largest: some subpixels may have the same gravitational value). Moreover, the gravitational values of N-class subpixels are sorted from smallest to largest: . In this research, the low attractiveness range for class M is [0, A M u ], and the low attractiveness range for class N is [0, A N v ]. Then, the swapping process can be expressed by the following formula: where A p i,j is the attractiveness of subpixel p i,j calculated by either Equation (1) or (3). Calculate the new total objective function values A P a,b and ∆E = A P a,b − A P a,b . If ∆E > 0, the new positions of the subpixels will be stored. If ∆E ≤ 0, compute the acceptance probability ε = exp(A P a,b − A P a,b )/T q at the current MSA T q . If ε > r, where r = rand (0,1), the new position of the subpixels will still be accepted. Otherwise, the locations of the subpixels will return to the previous state.
Step 3. Cooling system Once step 2 is accomplished, i = i + 1. If i = iter, the system will be cooled by the cooling function: Repeat steps 2 and 3. The optimization of one mixed pixel will be stopped when T q+1 < T stop . All mixed pixels will be optimized via steps 1 to 3.
Step 4. Randomly optimize the mixed pixels a second time The traditional PSA sequentially optimizes the mixed pixels, as shown in Figure 4a. However, the experimental results show that sequential optimization may limit the accuracy of SPM, which may be due to repeated one-directional optimization. To change the optimization direction, we randomly adjusted the order of the mixed pixels, as shown in Figure 4b and optimized the mixed pixels again following the same procedures described in steps 2 and 3. The experimental results in Section 3.4 show that randomly optimizing the mixed pixels again is better than sequentially optimizing them.  Step 5. Winner-takes-all strategy Since the SA algorithm uses binary optimization, only one type of land cover is optimized at a time. In each iteration, only the subpixels of one land cover type are set to 1 (pi,j = 1), and the pi,j for all the other subpixels are set to 0. To optimize all the land covers, steps 1 to 4 are repeated. After all land covers are optimized with the initial input proportional image, the corresponding class-specific attractiveness of the subpixels is calculated. Thus, each subpixel stored the attractiveness values corresponding to different types of land covers. Comparing the attractiveness values of all kinds of land covers within a subpixel, the subpixel will be classified according the land cover that has the largest attractiveness value 15.

Experiments and Analysis
To evaluate the proposed PSA_MSA algorithm, four different types of experiments were performed. Considering that SPM is only an image processing technology for proportion images, all tested images were classified and processed from high-resolution images. A mean filter is used to degrade these high-resolution images into low-resolution images with mixed pixels. The size of the mean filter is set to S × S, where S is the reconstruction scale. Thus, other errors can be avoided when generating proportion images, so that the reconstruction results are reliable and the positioning accuracy of different SPM methods are true and authentic. In the first three experiments, the PSA, the PSA_SA, and the DSAM were selected for comparison. In the fourth experiment, the influence of errors on the process of proportion acquisition in the SPM_MSA algorithm is tested. In the last experiment, the impact of the random optimization procedure of the PSA_MSA algorithm is discussed.
The overall accuracy is evaluated in terms of the percentage of correctly classified subpixels (PCC) and the classic Kappa coefficient [39]. The larger the PCC and Kappa values are, the more subpixels that are correctly classified and the better the reconstruction result. According to Merterns, the improved PCC' and Kappa' "are identical to PCC and Kappa except that they are calculated only for mixed pixels" [40]. The subpixels in pure pixels all belong to the same class and will undoubtedly increase both values of PCC and Kappa. Thus, the improved PCC' and Kappa' can better reflect the performance of different SPM methods than the original PCC and Kappa, respectively. Therefore, we also use the improved PCC' and Kappa' to evaluate the performance.

Artificial Shapes
In the first experiment, three different artificial shapes, as shown in Figure 5a, were used to test the algorithms. The reconstruction scale S was set to 10. In the SPM method based on SA, the parameters are set as follows: Tstart = 10 × S = 10 × 10, iter = 5, a = 0.8, and Tstop = 0.01. In the PSA_MSA model, u and v, which determine the range of subminimum attractiveness, were both set to 2. The results of the different algorithms are shown in Figure 5b-f. The corresponding evaluation criteria of the different methods are summarized in Table 1. It is obvious that the reconstruction results from directly using a conventional SA swapping algorithm are not as good as expected. This finding shows that the completely random character of the conventional SA algorithm does not fit the SPM problem very well. Comparing the results of the PSA and those obtained by the modified SA algorithms, the Step 5. Winner-takes-all strategy Since the SA algorithm uses binary optimization, only one type of land cover is optimized at a time. In each iteration, only the subpixels of one land cover type are set to 1 (p i,j = 1), and the p i,j for all the other subpixels are set to 0. To optimize all the land covers, steps 1 to 4 are repeated. After all land covers are optimized with the initial input proportional image, the corresponding class-specific attractiveness of the subpixels is calculated. Thus, each subpixel stored the attractiveness values corresponding to different types of land covers. Comparing the attractiveness values of all kinds of land covers within a subpixel, the subpixel will be classified according the land cover that has the largest attractiveness value 15.

Experiments and Analysis
To evaluate the proposed PSA_MSA algorithm, four different types of experiments were performed. Considering that SPM is only an image processing technology for proportion images, all tested images were classified and processed from high-resolution images. A mean filter is used to degrade these high-resolution images into low-resolution images with mixed pixels. The size of the mean filter is set to S × S, where S is the reconstruction scale. Thus, other errors can be avoided when generating proportion images, so that the reconstruction results are reliable and the positioning accuracy of different SPM methods are true and authentic. In the first three experiments, the PSA, the PSA_SA, and the DSAM were selected for comparison. In the fourth experiment, the influence of errors on the process of proportion acquisition in the SPM_MSA algorithm is tested. In the last experiment, the impact of the random optimization procedure of the PSA_MSA algorithm is discussed.
The overall accuracy is evaluated in terms of the percentage of correctly classified subpixels (PCC) and the classic Kappa coefficient [39]. The larger the PCC and Kappa values are, the more subpixels that are correctly classified and the better the reconstruction result. According to Merterns, the improved PCC and Kappa "are identical to PCC and Kappa except that they are calculated only for mixed pixels" [40]. The subpixels in pure pixels all belong to the same class and will undoubtedly increase both values of PCC and Kappa. Thus, the improved PCC and Kappa can better reflect the performance of different SPM methods than the original PCC and Kappa, respectively. Therefore, we also use the improved PCC and Kappa to evaluate the performance.

Artificial Shapes
In the first experiment, three different artificial shapes, as shown in Figure 5a, were used to test the algorithms. The reconstruction scale S was set to 10. In the SPM method based on SA, the parameters are set as follows: T start = 10 × S = 10 × 10, iter = 5, a = 0.8, and T stop = 0.01. In the PSA_MSA model, u and v, which determine the range of subminimum attractiveness, were both set to 2. The results of the different algorithms are shown in Figure 5b-f. The corresponding evaluation criteria of the different methods are summarized in Table 1. It is obvious that the reconstruction results from directly using a conventional SA swapping algorithm are not as good as expected. This finding shows that the completely random character of the conventional SA algorithm does not fit the SPM problem very well. Comparing the results of the PSA and those obtained by the modified SA algorithms, the modified SA algorithms obtain much better results, regardless of the type of shape in the image: circle, star, or letters. The PSA_MSA1 and PSA_MSA2 algorithms are based on our MSA algorithm with objective functions given by Equation (2) and Equation (4), respectively. The performances of these two algorithms are relatively close. The resulting images with letters show that the objective function based on Equation (4) is more suitable for more sophisticated mixed pixels with finer geometrical details and shapes than is the objective function based on Equation (2).
Sensors 2020, 20,1503 8 of 20 modified SA algorithms obtain much better results, regardless of the type of shape in the image: circle, star, or letters. The PSA_MSA1 and PSA_MSA2 algorithms are based on our MSA algorithm with objective functions given by Equation (2) and Equation (4), respectively. The performances of these two algorithms are relatively close. The resulting images with letters show that the objective function based on Equation (4) is more suitable for more sophisticated mixed pixels with finer geometrical details and shapes than is the objective function based on Equation (2).

Artificial Land Images
To demonstrate the spatial positioning ability of our algorithm on SPM, we applied the above algorithms to an artificial image containing four land cover types. The 450 × 450 artificial land image is shown in Figure 6a. Different shapes and boundaries are included in the artificial land image. The reconstruction scales were set to 8, 10, and 12. The parameters of the SA-based algorithms are the same as in Section 3.1. Figure 7a-e show the mapping results of the PSA, the PSA_SA algorithm, the DSAM algorithm, the PSA_MSA1 algorithm, and the PSA_MSA2 models at S = 12. It is obvious that the PSA_MSA1 and the PSA_MSA2 models reflect the spatial resolution of different classes better than the other three models did.

Artificial Land Images
To demonstrate the spatial positioning ability of our algorithm on SPM, we applied the above algorithms to an artificial image containing four land cover types. The 450 × 450 artificial land image is shown in Figure 6a. Different shapes and boundaries are included in the artificial land image. The reconstruction scales were set to 8, 10, and 12. The parameters of the SA-based algorithms are the same as in Section 3.1. Figure 7a-e show the mapping results of the PSA, the PSA_SA algorithm, the DSAM algorithm, the PSA_MSA1 algorithm, and the PSA_MSA2 models at S = 12. It is obvious that the PSA_MSA1 and the PSA_MSA2 models reflect the spatial resolution of different classes better than the other three models did.  The PCC, Kappa, improved PCC', and Kappa' values for different algorithms at the three reconstruction scales are summarized in Table 2. The larger the reconstruction scale is, the lower the accuracies of these methods because the mixed geometries are more complicated as the reconstruction scale increases; consequently, the uncertainty of the positioning process is inevitable.  Sensors 2020, 20, 1503 9 of 20

Artificial Land Images
To demonstrate the spatial positioning ability of our algorithm on SPM, we applied the above algorithms to an artificial image containing four land cover types. The 450 × 450 artificial land image is shown in Figure 6a. Different shapes and boundaries are included in the artificial land image. The reconstruction scales were set to 8, 10, and 12. The parameters of the SA-based algorithms are the same as in Section 3.1. Figure 7a-e show the mapping results of the PSA, the PSA_SA algorithm, the DSAM algorithm, the PSA_MSA1 algorithm, and the PSA_MSA2 models at S = 12. It is obvious that the PSA_MSA1 and the PSA_MSA2 models reflect the spatial resolution of different classes better than the other three models did.   Table 2. The larger the reconstruction scale is, the lower the accuracies of these methods because the mixed geometries are more complicated as the reconstruction scale increases; consequently, the uncertainty of the positioning process is inevitable.  The PCC, Kappa, improved PCC , and Kappa values for different algorithms at the three reconstruction scales are summarized in Table 2. The larger the reconstruction scale is, the lower the accuracies of these methods because the mixed geometries are more complicated as the reconstruction scale increases; consequently, the uncertainty of the positioning process is inevitable. At the same reconstruction scale, the PCC, Kappa, and improved PCC , Kappa values from the PSA, the PSA_SA and DSAM algorithm were lower than the values from the PSA_MSA1 and the PSA_MSA2 algorithms, which illustrates that the PSA_MSA algorithm performed better than the other three algorithms. Taking S = 12 as an example, the improved Kappa value from the PSA_MSA1 and the PSA_MSA2 algorithms reached 94.77 and 94.99, respectively. The improved Kappa value from the PSA_MSA1 was 12.79%, which is 20.80% higher than that of the PSA and the PSA_SA algorithm. Additionally, it is 4.86% higher than that of DSAM as well. Meanwhile, the improved Kappa value from the PSA_MSA2 algorithm was 13.06%, which is 21.08% higher than that of the PSA and the PSA_SA algorithm. In addition, it is 5.10% higher than that of DSAM as well. Comparing the results of the PSA_MSA1 and the PSA_MSA2 algorithms, the PSA_MSA1 algorithm performs slightly better than the PSA_MSA2 algorithm at low reconstruction scales. As the reconstruction scale increases, the PSA_MSA2 algorithm performs slightly better than the PSA_MSA1 algorithm. This finding further confirms that Equation (4) is more suitable for low-resolution images with more sophisticated geometries than Equation (2).
To further show that the proposed PSA_MSA algorithm is superior to PSA, PSA_SA, and DSAM, we generated a synthetic data set which contains 108 artificial land images of size 240 × 240 and four kinds of classes, including various shapes. The reconstruction scales were set to 8, 10, and 12. Using PCC and Kappa as evaluation indexes, the results are plotted as curves in Figures 8 and 9. Equations (6) and (7) are used to compare the accuracies of PSA_MSA1 and PSA_MSA2 with those of PSA, PSA_SA, and DSAM, respectively. In this experiment, num = 108.
Sensors 2020, 20, 1503 10 of 20 At the same reconstruction scale, the PCC, Kappa, and improved PCC', Kappa' values from the PSA, the PSA_SA and DSAM algorithm were lower than the values from the PSA_MSA1 and the PSA_MSA2 algorithms, which illustrates that the PSA_MSA algorithm performed better than the other three algorithms. Taking S = 12 as an example, the improved Kappa' value from the PSA_MSA1 and the PSA_MSA2 algorithms reached 94.77 and 94.99, respectively. The improved Kappa' value from the PSA_MSA1 was 12.79%, which is 20.80% higher than that of the PSA and the PSA_SA algorithm. Additionally, it is 4.86% higher than that of DSAM as well. Meanwhile, the improved Kappa' value from the PSA_MSA2 algorithm was 13.06%, which is 21.08% higher than that of the PSA and the PSA_SA algorithm. In addition, it is 5.10% higher than that of DSAM as well. Comparing the results of the PSA_MSA1 and the PSA_MSA2 algorithms, the PSA_MSA1 algorithm performs slightly better than the PSA_MSA2 algorithm at low reconstruction scales. As the reconstruction scale increases, the PSA_MSA2 algorithm performs slightly better than the PSA_MSA1 algorithm. This finding further confirms that Equation (4) is more suitable for low-resolution images with more sophisticated geometries than Equation (2).
To further show that the proposed PSA_MSA algorithm is superior to PSA, PSA_SA, and DSAM, we generated a synthetic data set which contains 108 artificial land images of size 240 × 240 and four kinds of classes, including various shapes. The reconstruction scales were set to 8, 10, and 12. Using PCC' and Kappa' as evaluation indexes, the results are plotted as curves in Figures 8 and 9. Equations (6) and (7) are used to compare the accuracies of PSA_MSA1 and PSA_MSA2 with those of PSA, PSA_SA, and DSAM, respectively. In this experiment, num = 108.   As shown in Tables 3 and 4, the larger the reconstruction scale is, the better the performance of PSA_MSA compared to those of PSA, PSA_SA, and DSAM. The results further confirmed that PSA_MSA achieves a better accuracy than PSA, PSA_SA, and DSAM in SPM.  As shown in Tables 3 and 4, the larger the reconstruction scale is, the better the performance of PSA_MSA compared to those of PSA, PSA_SA, and DSAM. The results further confirmed that PSA_MSA achieves a better accuracy than PSA, PSA_SA, and DSAM in SPM.

Real land Image
In this section, to further demonstrate the advantages of the PSA_MSA algorithm, two real land images are selected as the experimental objects. One is a highly developed area with high spatial autocorrelation and is marked as real land image1, while the other, which is marked as real land image2, is a region with low spatial autocorrelation that has not been highly developed.
Real land image1 is a remote sensing image of Nanjing, Jiangsu Province. The image was taken by the Landsat 8 satellite, which was launched in 2013. The 30-m reference image was taken on October 9, 2017, and contains a total of 396 × 396 pixels. After the image was processed by using K-means classification [33,41], three different classes which shown in Figure 10a were obtained: water, vegetation and soil. Real land image2 is a remote sensing image of Ningbo, Zhejiang Province. The 30-m reference image was taken by the Landsat 8 satellite on August 24, 2017, and contains a total of 398 × 398 pixels. The data can be downloaded from the following website: http://www.gscloud.cn/. The S × S mean filter was used to degrade these two high-resolution images. The reconstruction scale was set to 8 and 10. The parameters of the SA-based algorithms are mentioned in Section 3.1.
October 9, 2017, and contains a total of 396 × 396 pixels. After the image was processed by using Kmeans classification [33,41], three different classes which shown in Figure 10(a) were obtained: water, vegetation and soil. Real land image2 is a remote sensing image of Ningbo, Zhejiang Province. The 30-m reference image was taken by the Landsat 8 satellite on August 24, 2017, and contains a total of 398 × 398 pixels. We also used the K-means method to obtain three different classes in this image Figure 12(a): water, vegetation, and soil. The data can be downloaded from the following website: http://www.gscloud.cn/. The S × S mean filter was used to degrade these two high-resolution images. The reconstruction scale was set to 8 and 10. The parameters of the SA-based algorithms are mentioned in Section 3.1.
(1) Real land image1 Figure 11a-f display the SPM results of the PSA, the PSA_SA, the DSAM algorithm, the PSA_MSA1 algorithm, and the PSA_MSA2 algorithm at S = 8, respectively. Although there were some jagged boundaries for different land covers, the SPM results of the MSA-based algorithms were better than the results of other algorithms. To further demonstrate the superiority of the MSA-based algorithms, a local area was selected and enlarged to show the mapping results. Obviously, the PSA_MSA algorithms reconstructed the target image much more effectively than the other methods.  Table 5 summarizes the SPM results of various models. Both the PCC' and Kappa' values from the PSA_MSA1 and PSA_MSA2 algorithms are much better than those values from the other algorithms. Taking S = 10 as an example, the improved Kappa' value from the PSA_MSA1 and PSA_MSA2 algorithms reached 71.42 and 71.44, respectively. The improved Kappa' value from the PSA_MSA1 algorithm was 11.28%, which is 18.58% higher than that of the PSA and the PSA_SA algorithm. Moreover, it is 4.54% higher than that of DSAM as well. Additionally, the improved Kappa' value from the PSA_MSA2 algorithm was 11.31%, which is 18.61% higher than that of the PSA and the PSA_SA algorithm. It is 4.57% higher than that of DSAM as well. The PSA_MSA2 algorithm performed slightly better than the PSA_MSA1 algorithm in the case of complex mixing land cover pixels.
The SPM accuracy (improved PCC') for each land cover type from the five methods is also listed in Table 6. All land covers achieved the highest accuracy levels with the PSA_MSA algorithms. (1) Real land image1 Figure 11a-f display the SPM results of the PSA, the PSA_SA, the DSAM algorithm, the PSA_MSA1 algorithm, and the PSA_MSA2 algorithm at S = 8, respectively. Although there were some jagged boundaries for different land covers, the SPM results of the MSA-based algorithms were better than the results of other algorithms. To further demonstrate the superiority of the MSA-based algorithms, a local area was selected and enlarged to show the mapping results. Obviously, the PSA_MSA algorithms reconstructed the target image much more effectively than the other methods.     11.28%, which is 18.58% higher than that of the PSA and the PSA_SA algorithm. Moreover, it is 4.54% higher than that of DSAM as well. Additionally, the improved Kappa value from the PSA_MSA2 algorithm was 11.31%, which is 18.61% higher than that of the PSA and the PSA_SA algorithm. It is 4.57% higher than that of DSAM as well. The PSA_MSA2 algorithm performed slightly better than the PSA_MSA1 algorithm in the case of complex mixing land cover pixels. The SPM accuracy (improved PCC ) for each land cover type from the five methods is also listed in Table 6. All land covers achieved the highest accuracy levels with the PSA_MSA algorithms. (2) Real land image2 We also used the K-means method to obtain three different classes in this image Figure 12a: water, vegetation, and soil. The SPM results of the PSA, the PSA_SA, the DSAM algorithm, the PSA_MSA1 algorithm and the PSA_MSA2 algorithm at S = 8 are displayed in Figure 13a-e, respectively. The accuracies of all the algorithms for real land image2 are lower than the result of SPM for real land image1. This shows that these algorithms are more suitable for images with high spatial autocorrelation.  As seen in Table 7, if the land covers in the image are trivial and lack spatial correlation, then the mapping accuracy of all algorithms based on spatial dependence theory will decrease. Although the accuracies of all the algorithms decreased, it can be seen from Table 7 that both the improved PCC' and Kappa' values from the PSA_MSA1 and PSA_MSA2 algorithms are still better than the values from the other algorithms. However, it should be pointed out that if the spatial distribution of the image is too scattered, then PSA_MSA will also fail. Thus, for the image with weak spatial correlation, it is necessary to improve the accuracy of proportional image acquisition and introduce variation function to subpixel mapping. The SPM accuracies (improved PCC') for each land cover type from the five methods are also listed in Table 8. Vegetation achieved the highest accuracies with the PSA_MSA algorithm.  As seen in Table 7, if the land covers in the image are trivial and lack spatial correlation, then the mapping accuracy of all algorithms based on spatial dependence theory will decrease. Although the accuracies of all the algorithms decreased, it can be seen from Table 7 that both the improved PCC' and Kappa' values from the PSA_MSA1 and PSA_MSA2 algorithms are still better than the values from the other algorithms. However, it should be pointed out that if the spatial distribution of the image is too scattered, then PSA_MSA will also fail. Thus, for the image with weak spatial correlation, it is necessary to improve the accuracy of proportional image acquisition and introduce variation function to subpixel mapping. The SPM accuracies (improved PCC') for each land cover type from the five methods are also listed in Table 8. Vegetation achieved the highest accuracies with the PSA_MSA algorithm. As seen in Table 7, if the land covers in the image are trivial and lack spatial correlation, then the mapping accuracy of all algorithms based on spatial dependence theory will decrease. Although the accuracies of all the algorithms decreased, it can be seen from Table 7 that both the improved PCC and Kappa values from the PSA_MSA1 and PSA_MSA2 algorithms are still better than the values from the other algorithms. However, it should be pointed out that if the spatial distribution of the image is too scattered, then PSA_MSA will also fail. Thus, for the image with weak spatial correlation, it is necessary to improve the accuracy of proportional image acquisition and introduce variation function to subpixel mapping. The SPM accuracies (improved PCC ) for each land cover type from the five methods are also listed in Table 8. Vegetation achieved the highest accuracies with the PSA_MSA algorithm.

Influence of Errors on the Process of Proportion Acquisition in the SPM_MSA Algorithm
It is unavoidable that there will be some errors in the process of obtaining end element proportions caused by spectral unmixing. To further investigate the influence of proportion errors on the PSA_MSA algorithm, we introduce a random error in the process of generating the proportional image, which creates a certain number of errors between the proportional image applied to the PSA_MSA algorithm and the actual proportional image. Real land image1 was selected as the test object. We added 5% and 10% random errors in the proportion image. Figure 14 shows proportion images with an increase of 10% error and containing no pure pixels. The SPM result can be seen in Figure 15. Since the proportion errors are added to each mixed pixel, PCC is equivalent to the improved PCC and Kappa is equivalent to the improved Kappa . Thus, PCC and Kappa are selected as the evaluation indexes for comparison with the case of 0 error. Tables 9 and 10 summarize the subpixel mapping accuracies of the different algorithms under different errors. It is not difficult to see that the subpixel mapping accuracies of all the algorithms decrease with increasing errors. This shows that errors in the input data will reduce the SPM accuracy. As shown in Figure 15, all the algorithms maximize the spatial correlations of the proportion errors of different land covers existing in the mixed pixels, thereby aggregating them and badly affecting the accuracy of the final subpixel mapping. Tables 9 and 10 show that the lower the error of the proportion value is, the higher the accuracy of the PSA_MSA algorithm is. This shows the effectiveness of the proposed PSA_MSA algorithm. In other words, if the accuracy of proportion data can be improved, the advantage of PSA_MSA will be more obvious. Compared with PSA and PSA_SA, both DSAM and PSA_MSA can provide SPM results with the best visual and numerical effects.

Influence of Errors on the Process of Proportion Acquisition in the SPM_MSA Algorithm
It is unavoidable that there will be some errors in the process of obtaining end element proportions caused by spectral unmixing. To further investigate the influence of proportion errors on the PSA_MSA algorithm, we introduce a random error in the process of generating the proportional image, which creates a certain number of errors between the proportional image applied to the PSA_MSA algorithm and the actual proportional image. Real land image1 was selected as the test object. We added 5% and 10% random errors in the proportion image. Figure 14 shows proportion images with an increase of 10% error and containing no pure pixels. The SPM result can be seen in Figure 15. Since the proportion errors are added to each mixed pixel, PCC is equivalent to the improved PCC' and Kappa is equivalent to the improved Kappa'. Thus, PCC and Kappa are selected as the evaluation indexes for comparison with the case of 0 error. Tables 9 and 10 summarize the subpixel mapping accuracies of the different algorithms under different errors. It is not difficult to see that the subpixel mapping accuracies of all the algorithms decrease with increasing errors. This shows that errors in the input data will reduce the SPM accuracy. As shown in Figure 15, all the algorithms maximize the spatial correlations of the proportion errors of different land covers existing in the mixed pixels, thereby aggregating them and badly affecting the accuracy of the final subpixel mapping. Tables 9 and 10 show that the lower the error of the proportion value is, the higher the accuracy of the PSA_MSA algorithm is. This shows the effectiveness of the proposed PSA_MSA algorithm. In other words, if the accuracy of proportion data can be improved, the advantage of PSA_MSA will be more obvious. Compared with PSA and PSA_SA, both DSAM and PSA_MSA can provide SPM results with the best visual and numerical effects.

Stochastic Optimization of Mixed Pixels in the PSA_MSA Algorithm
In step 4 of the proposed PSA_MSA algorithm, we optimized the mixed pixels in a random sequence, which is called stochastic optimization. In traditional optimization, the mixed pixels are sequentially processed by rows or columns, as shown in Figure 4a. Stochastic optimization can change the optimization direction. To analyze the effect stochastic optimization for the PSA_MSA algorithm, the two optimization methods were implemented in the secondary optimization phase of the PSA_MSA algorithm, and the positioning accuracy results were compared. The iteration times for the two optimization methods were the same, and the total optimizations accepted by a single mixed pixel were also the same. Circles and stars were selected as the test objects. Table 11 shows the results after setting the reconstruction scale to 8, and Table 12 shows the results after setting the reconstruction scale to 10. Stochastic optimization can achieve a better SPM accuracy than sequential optimization.

Stochastic Optimization of Mixed Pixels in the PSA_MSA Algorithm
In step 4 of the proposed PSA_MSA algorithm, we optimized the mixed pixels in a random sequence, which is called stochastic optimization. In traditional optimization, the mixed pixels are sequentially processed by rows or columns, as shown in Figure 4a. Stochastic optimization can change the optimization direction. To analyze the effect stochastic optimization for the PSA_MSA algorithm, the two optimization methods were implemented in the secondary optimization phase of the PSA_MSA algorithm, and the positioning accuracy results were compared. The iteration times for the two optimization methods were the same, and the total optimizations accepted by a single mixed pixel were also the same. Circles and stars were selected as the test objects. Table 11 shows the results after setting the reconstruction scale to 8, and Table 12 shows the results after setting the reconstruction scale to 10. Stochastic optimization can achieve a better SPM accuracy than sequential optimization.

Discussion
The primary different between the PSA_MSA method with the PSA_SA method is that the selection range for the subpixels to be swapped is narrower in the PSA_MSA method. The results shown in Section 3 demonstrate that the SPM positioning accuracy is significantly improved. In Section 3,

Discussion
The primary different between the PSA_MSA method with the PSA_SA method is that the selection range for the subpixels to be swapped is narrower in the PSA_MSA method. The results shown in Section 3 demonstrate that the SPM positioning accuracy is significantly improved. In Section 3, we randomly selected an M-class subpixel of in the gravitational range of [0, A M 2 ] and an N-class subpixel in the gravitational range of [0, A N

Discussion
The primary different between the PSA_MSA method with the PSA_SA method is that the selection range for the subpixels to be swapped is narrower in the PSA_MSA method. The results shown in Section 3 demonstrate that the SPM positioning accuracy is significantly improved. In Section 3, we randomly selected an M-class subpixel of in the gravitational range of [0, A M 2 ] and an N-class subpixel in the gravitational range of [0, A N 2 ] for swapping. In fact, other low attractiveness ranges can be selected, such as [0, A M 3 ] and [0, A N 3 ] or [0, A M m ] and [0, A N n ]. To analyze the effect of the low attractiveness range, six different low attractiveness ranges were selected for the PSA_MSA algorithm: from [0, A M 1 ] and [0, A N 1 ] to [0, A M 6 ] and [0, A N 6 ], respectively. Therefore, in the experiments, m = n = w (w= 1, 2, …,6). Circles and stars were selected as the test objects, and the SPM accuracies are shown in Figures 16 and 17.
The above results demonstrate that the SPM accuracy is better when w = 2 or 3, which corresponds to [0, A M 2 ] and [0, A N 2 ] and [0, A M 3 ] and [0, A N 3 ], respectively. As the range of low attractiveness widens, the SPM accuracy decreases. However, the optimal range for subpixel swapping might be different in different scenarios. In this research, all the experiments in Section 3 are performed with w = 2.   The above results demonstrate that the SPM accuracy is better when w = 2 or 3, which corresponds to [0, A M 2 ] and [0, A N 2 ] and [0, A M 3 ] and [0, A N 3 ], respectively. As the range of low attractiveness widens, the SPM accuracy decreases. However, the optimal range for subpixel swapping might be different in different scenarios. In this research, all the experiments in Section 3 are performed with w = 2.

Conclusions
SPM methods can improve the spatial resolution of spectral images. The traditional PSA model is prone to fall into a local optimum during subpixel swapping. To solve this problem, we proposed a PSA_MSA method that uses SA to randomly swap two subpixels from two different classes. These two subpixels are selected from subpixels with class-specific attractiveness values in certain low ranges. This SA algorithm enables subpixel swapping to escape the local optimum solutions. Only the subpixels with low class-specific attractiveness values may be selected for swapping. Two spatial correlation models were applied to calculate the attractiveness coefficients and compare the mapping results. Experiments were performed with different shapes and land cover images. Compared with other models, such as the PSA, the PSA_SA, and the DSAM algorithm, the proposed PSA_MSA performed better than the others did in the reconstruction of spectral images. The improved PCC and Kappa values also showed that the proposed PSA_MSA algorithm can generate higher accuracy SPM results than the other algorithms. The PSA_MSA method involved random selection. The experiments also confirmed that a randomized optimization sequence can also improve the mapping accuracy.
The effect of errors in the input proportion images were discussed. The results showed that the larger values of proportion errors lead to low mapping accuracy of SPM algorithms. These proportion errors are mistaken as wrong end-members by the algorithms. SPM algorithms will transform the wrong proportion information of the end-members into spatial information. To reduce the effects of the proportion errors, it is very important to obtain more accurate proportion information in the spectral unmixing procedure. In addition, the PSA_MSA algorithm needs a large number of iterations to obtain global solution. It is meaningful to look for an improved method to shorten the running time of PSA_MSA algorithm. We are considering the following work in the future: First, we aim to study the spectral unmixing algorithms and improve the estimation accuracy to reduce the errors in the proportion images. Second, we aim to further optimize the mapping result, the spectral unmixing procedure, and end-member mapping process while combining all three in an integrated framework. Third, further analysis will be conducted regarding factors that may affect the performance of subpixel mapping algorithms based on spatial dependence theory, such as the spatial resolutions of remote sensed images, the number of end-members in a mixed pixel, and the scattering distribution of the end-members.
Author Contributions: L.S. proposed the concept of the algorithm and wrote the paper; Y.X. designed the algorithm, Y.X. and J.Y. performed the experiments and analyzed the results; Y.Y. adjusted the algorithm and edited the framework of the experiment and the paper. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.