A Novel Unsupervised Classiﬁcation Method for Sandy Land Using Fully Polarimetric SAR Data

: Classiﬁcation based on polarimetric synthetic aperture radar (PolSAR) images is an emerg-ing technology, and recent years have seen the introduction of various classiﬁcation methods that have been proven to be effective to identify typical features of many terrain types. Among the many regions of the study, the Hunshandake Sandy Land in Inner Mongolia, China stands out for its vast area of sandy land, variety of ground objects, and intricate structure, with more irregular characteristics than conventional land cover. Accounting for the particular surface features of the Hunshandake Sandy Land, an unsupervised classiﬁcation method based on new decomposition and large-scale spectral clustering with superpixels (ND-LSC) is proposed in this study. Firstly, the polarization scattering parameters are extracted through a new decomposition, rather than other decomposition approaches, which gives rise to more accurate feature vector estimate. Secondly, a large-scale spectral clustering is applied as appropriate to meet the massive land and complex terrain. More speciﬁcally, this involves a beginning sub-step of superpixels generation via the Adaptive Simple Linear Iterative Clustering (ASLIC) algorithm when the feature vector combined with the spatial coordinate information are employed as input, and subsequently a sub-step of representative points selection as well as bipartite graph formation, followed by the spectral clustering algorithm to complete the classiﬁcation task. Finally, testing and analysis are conducted on the RADARSAT-2 fully PolSAR dataset acquired over the Hunshandake Sandy Land in 2016. Both qualitative and quantitative experiments compared with several classiﬁcation methods are conducted to show that proposed method can signiﬁcantly improve performance on classiﬁcation. clustering with superpixels. We propose improvements from these two aspects to classify the sandy land.


Introduction
In recent years, desertification of grasslands has become increasingly severe due to continuous drought and grassland overload in Inner Mongolia Autonomous Region, China. According to research, Hunshandake Sandy Land, has become one of the main dust sources, so it is of great significance to monitor sandy land [1]. The development of spaceborne remote sensing technology and a massive increase in monitoring data provide many effective means to monitor the sandy land echo-environment [1,2]. The all-time all-weather features of active microwave synthetic aperture radar (SAR) are a particular advantage over optical sensors, and SAR can penetrate ground and vegetation with multi-polar observations [3,4]. The polarimetric SAR (PolSAR) which can obtain relatively complete polarimetric scattering information through four different channels (HH, HV, VH and VV) is referred to as a fully PolSAR [5]. Theory and practice show that the polarization feature of electromagnetic echoes is sensitive to the shape, texture, and other physical features of targets, and the differences in polarization signatures help natural images, which has been successfully applied for PolSAR images [25,29,33,34]. The algorithm can generate superpixels with compact regions, relatively uniform sizes, and maintain the region boundaries. In this paper, the PolSAR image is segmented by Adaptive Simple Linear Iteration Clustering (ASLIC), and the image is divided into different region blocks. These area blocks are used as the input sample points of the spectrum clustering algorithm. Compared with directly using the original image, the amount of data processing is greatly reduced, the problems of insufficient memory and slow processing speed are improved [35]. Furthermore, as for clustering algorithms, such as k-means and expectation maximization (EM), must assume a convex spherical sample space, and they tend to fall into local optima [23][24][25][26]. Fortunately, the spectral clustering algorithm, which is based on graph theory, can transform clustering into a graph partitioning problem. It is carried out by k-means clustering in a projection space, in which the transformation matrix is calculated by eigenvectors of Gaussian similarity matrix of samples [24]. Such methods are attracting considerable attention in the field of polarized SAR processing.
In this study, Hunshandake Sandy Land stands out for its vast area of sandy land, a great variety of ground objects, and a complex, intricate structure with irregular characteristics compared with conventional land cover types. In order to mitigate the problem and improve final classification accuracy, starting from the two key steps of unsupervised classification, an unsupervised classification method based on new decomposition and large-scale spectral clustering with superpixels (ND-LSC) is proposed. Figure 1 depicts the proposed architecture. There are two essential parts for the classification method. The first part is the extraction of polarized scattering power by ND which contribute significantly to understand the polarimetric scattering mechanisms of sandy land. The second part is the large-scale spectral clustering (LSC), which is a classification algorithm based on regional level. The algorithm, i.e., scattering power as input vs superpixels as basic unit, is tested on the sandy land. Specifically, we utilize the advantages of new decomposition to extract polarization scattering power (m s , m d , and m v ), which is used to construct the feature vector X = [m s , m d , m v ], and then superpixels are generated by the ASLIC algorithm when the feature vector combined with the spatial coordinate information are employed as input, followed by spectral clustering to complete the classification. There are several advantages of the proposed method. The main novel contributions of this work are: Firstly, it is a classification structure for areas with vast area, variety of ground objects and intricate structure.
Secondly, these two aspects are crucial to complete the final classification through spectral clustering. On the one hand, the new decomposition method helps to interpret the polarization scattering characteristics of ground objects. On the other hand, the classification algorithm takes the superpixel as the unit. In LSC, it can make full use of the polarization information and spatial information of the target to generate superpixels. These two aspects are crucial to complete the final classification through spectral clustering.
The structure of this paper is organized as follows. In Section 2, the descriptions of the Hunshadake Sandy Land and parameters of RADARSAT-2 fully PolSAR were presented. Then, the details of the proposed ND-LSC unsupervised classification are introduced in Section 3, and the flowchart of the method is given at the end of this section. In Section 4, on the one hand, the classification performance of ND-LSC and other classification methods are compared. On the other hand, it analyzed the influence of decomposition parameters and superpixels segmentation on classification result. Finally, discussion and conclusions are drawn in Sections 5 and 6, respectively.

Study Area
Hunshandake Sandy Land, is located at 42 • 32 -42 • 50 N, 115 • 49 -116 • 12 E in the south of Xilingol League, Inner Mongolia, China. It is a semi-temperate grassland with a semi-arid continental climate. The altitude ranges from 1150 to 1500 m rising slowly from northwest to southeast. The edge of the sandy land is hilly and denuded, and there are dunes, basins, lakes, and denuded plateaus in this area. Sand dunes can be divided into fixed, semi-fixed, and mobile dunes from the perspective of vegetation coverage. The vegetation of fixed dunes is mainly herbaceous, and the proportions of irrigation are relatively high. Semi-fixed sand dunes distributed among fixed sand dunes always present speckle shapes in SAR images. A mobile sand dune can move, as the effect of the wind for the surface vegetation is rare. There are five types of sandy land in the area. Fixed sand accounts for 29.23%, semi-fixed sand for 29.52%, mobile sand for 15.23%, semi-mobile sandy land for 23.27%, and saline-alkali land for only 2.55%. There is almost no Gobi sand in Hunshandake Sandy Land [35], as shown in Figure 2a. Polarimetric SAR data are obtained from the RADARSAT-2 satellite, with data parameters as shown in Table 1. The corresponding PolSAR images are shown in Figure 2b-d. According to surveys, three main types of phenomena occur or have occurred in Hunshandake Sandy Land: (1) soil and water are severely lost, (2) the degree of desertification increases, and (3) grass pasture is severely degraded. Therefore, it is urgent to constantly monitor the change of grassland features. The data used in this article were collected in September 2013, August 2016, and September 2017.  Hunshandake Sandy Land is famous for its unique characteristics, such as a large area, various types of features, and intricate structures. As shown in Figure 3, the details of seven typical features.

Methodology
An unsupervised classification method combined with new decomposition and largescale spectral clustering with superpixels is proposed for Hunshadake Sandy Land classification. The method consists of data preprocessing, new decomposition (ND), superpixel segmentation by adaptive simple linear iterative clustering (ASLIC), and large-scale spectral clustering (LSC) classification. The workflow of the proposed ND-LSC is shown in Figure 4. The polarization parameters (m s , m d , m v ) are more accurately extracted from the ND, which is used to construct the feature vector X = [m s , m d , m v ]. The application of LSC is more in line with the characteristics of its large area. The feature vector X and spatial coordinate information constitute the input vector of ASLIC, which is used to generate superpixels. Representative points are selected from superpixels, and the bipartite graph between the representative points is formed. According to the principle of the spectral clustering algorithm, the classification of PolSAR is completed. It can be seen from the flowchart that the classification method has two three parts, (1) extracting polarization parameters from decomposition (2) spectral clustering with superpixels. We propose improvements from these two aspects to classify the sandy land.

Data Preprocessing
The SAR pre-processing procedure is performed in three steps with SARscape software (PolSARpro 5.3.0). These include calibration, registering, and filtering. The PolSAR single look complex (SLC) image data can be geometrically calibrated and registered by the NEST tool. Then the Lee filter is used to suppress inherent speckle noise, and the corresponding polarimetric coherency matrix T 3 for every pixel of data is generated. Finally, the area of interest is selected.
As for fully PolSAR image, the acquired scattering matrix in H-V polarization basis can be written as [5]: where the different elements in matrix S represent the backscattering coefficients of the different polarimetric channels. S HH and S VV are the common-polarized echo channels, and S HV and S VH are the cross-polarized echo channels. According to the reciprocity principle, the two cross-polarizations are considered equal, that is, S HV = S VH , and the scattering target vector obtained from the matrix S can be expressed as follows: Then the corresponding coherency matrix is (3) [11]: T 3 contains nine independent parameters, which is relevant to the physical properties of the observed targets. Generally speaking, T 11 is relevant to surface scattering mechanism, T 22 is relevant to double-scattering mechanism, and T 33 is relevant to volume scattering mechanism. T 13 is relevant to cross-polarization power generated by the coupling between surface and volume scattering mechanisms. T 23 is relevant to cross-polarization power generated by the coupling between double-bounce and volume scattering mechanism. Elimination of the cross-polarization power cannot only optimize the PolSAR coherency matrix T 3 , but can help approach the reflection symmetry condition.

New Decomposition
In order to effectively obtain the polarization parameters of Hunshandake Sandy Land to improve the final classification results, a new decomposition is presented based on [7,17] in this paper. The decomposition eliminates cross-polarization at the bottom power to reduce the influence of the volume scattering overestimation problem, followed by the hybrid Freeman/Eigenvalue decomposition based on adaptive model to extract the polarization scattering power.

Removal
Orientation angle compensation (OAC) [15,16] and phase angle rotation (PAR) [18] are common polarimetric target decomposition data processing methods. They are used to remove T 23 and T 13 , which can increase the power of the common polarization term, decrease the power of the cross-polarization term, minimize the influence of different orientation angles on polarimetric decomposition results, Thus, the coherency matrix T 3 is rotated as follows: where *T denotes complex conjugation and transposition, and the rotation matrices R(ϕ 1 ) and R(ϕ 2 ) are given by [21]: After this rotation, the eliminated T 13 of coherency matrix T 3 can be eliminated by Equations (4) and (5) completely.
After OAC, the real part of T 13 becomes zero, and its imaginary part can be set to zero via PAR [21], that is, Similarly, T 23 represents a cross term of double-bounce and volume scattering, which can be completely eliminated by OAC and PAR, that is, The rotation matrices R(ϕ 1 ) and R(ϕ 2 ) are given by [18]. The choice of rotation transformation depends on the form of the dominant scattering mechanism. It can be determined according to the value of T 11 − T 22 , which contributes to removing the maximum of the cross-polarization power generated by the coupling between common polarization and the cross-polarization scattering mechanism. The coherent matrix after two unitary rotations is as follows:

Adaptive Model-Based Decomposition
After two unitary rotations, the coherency matrix T 3 ' can be decomposed as follows: where T s , T d and T v are scattering models of surface, double-bounce, and volume scattering, respectively, with corresponding power coefficients m s , m d and m v . The surface scattering model T s and double-bounce scattering model T d are represented by eigenvectors, and can be expressed as follows [13]: where α s and α d depend on the two dielectric constants of the surface, reflector and the angle of incidence, and α s < π/4, α d > π/4. ϕ s and ϕ d are the scattering phases for surface scattering and double-bounce scattering mechanisms, respectively. Under the conditions α s + α d = π/2 and ϕ d − ϕ s = ±π, these two components T s and T d have orthogonality with each other. This reduces the number of unknowns in (11), which can be rewritten as: The adaptive volume scattering model is the main distinction of improved hybrid decomposition, and can be described as follows: The object is judged to be located at artificial area when double-bounce scattering is dominant scattering mechanisms. Under the condition, we consider this physical situation in more detail and assume for more accurate modeling that volume scattering mainly comes from the dihedral structure of the artificial target. The model is expressed as follows [13]: The object is judged to be located in a natural area when surface scattering is dominant. At this time, volume scattering mainly comes from a cloud of randomly oriented dipole scatter such as vegetation, in which case we adopt the generalized volume scattering model [16]: where γ = |S HH | 2 / |S VV | 2 , and γ is the ratio of the horizontal and vertical polarization components.
To determine the scattering power coefficients of m v , m s and m d , first the volume scattering power m v coefficient is calculated as [17]: Then, Equation (14) can be written as follows: Finally, the remaining scattering power coefficient m s and m d can be calculated as the eigenvalues of remainder matrix T sd on the basis of dominant scattering mechanism (λ 1sd > λ 2sd ). In case of T 11 − T 22 > 0, m s and m d can be determined as: In the case of T 11 − T 22 < 0, m s and m d can be determined as:

LSC Classification
Spectral clustering algorithms usually involve two large steps (1) to construct a similarity matrix (2) to solve the corresponding eigenvalue problem on the Laplace. However, accounting for the particular surface features of the Hunshandake Sandy Land, the implementation of the algorithm requires a lot of memory and affects the classification accuracy. Therefore, on the basis of obtaining the polarization scattering parameters (m s , m d , m v ), this paper generates representative points through superpixels, and subsequently a sub-step of representative points selection as well as bipartite graph formation, followed by the spectral clustering algorithm to complete the classification task.

Superpixels Segmentation
The ASLIC algorithm is utilized to generate superpixels in this paper, and it is developed based on the SLIC algorithm. The basic principle of SLIC is to use k-means for segmentation in a local area. In the beginning, some seeds are uniformly selected on the image plane at intervals Ns as the initial clustering center of the superpixels. The formula for the distance between a pixel and its adjacent superpixel region in the PolSAR image is generally expressed as [32]: where d xy is the Euclidean distance between the image plane coordinate points i and j, and d W the Wishart distance. Since the Wishart distance may be negative, it is not a good distance measure. In this paper, Wishart distance will be replaced by Bartlett distance as: where T i and ∑ j are the respective sample coherence matrices of pixel S i and superpixel region S pi . Polarization parameters (m s , m d , m v ) can be extracted more accurately from ND, and are used to construct the vector X= [m s , m d , m v ] as the feature vector of the coherence matrix.
Since the values of d W and d XY changes with different images, and also with regions, D SLIC is modified to adaptively balance the distance between the pixel values and the distance in the image plane space, and Equation (22) can be modified as follows: The ranges of the two items in Equation (24) are very similar, and it is easy to determine the contribution of the two components. Consequently, ASLIC with the advantage of adaptive distance formula is selected to generate superpixel in this paper.

LSC Classification
The idea of spectral clustering algorithms is derived from the theory of spectral graph partitioning. If a dataset is regarded as a set of all vertices in a graph, and the weights of edges between vertices represent the similarities between data, then spectral clustering is a method to divide the graph by a weight. Then, to find the optimal solution of the graph partitioning criterion becomes the essential problem of spectral clustering. According to the Rayleigh-Ritz theory, several eigenvectors of similarity matrix W are optimal partitionapproximate solutions. That is, the original optimal segmentation problem is converted to the spectral decomposition of the similarity matrix or Laplacian matrix. Therefore, the key to classification for a PolSAR image based on spectral clustering is to obtain the similarity matrix W, expressed as follows [25]: where s i , s j are sample of data; d(s i , s j ) is the distance between superpixel, and usually is the Euclidean distance; and σ is a scale parameter. We construct a normalization Laplacian matrix L rw by similarity matrix W as [32]: where L = D − W, and D is the degree matrix, which consists of the corner elements of W. Based on the symmetric Wishart distance measure as Equation (23), the similarity matrix of spectral clustering can be defined as shown: where T pi and T pj are the coherency matrix of superpixels (p i , p j ), and d w is the symmetric Wishart distance. The similarity matrix of spectral clustering based on central features can be defined as shown: where F pi and F pj are the central features vector (X = [m s , m d , m v ]) of superpixels (p i , p j ), and d is the cosine distance. The similarity matrix Z for LSC is defined as follows: The procedure of the LSC in this paper can then be described as follows: 1. Segment the PolSAR image into superpixels with ASLIC.

2.
Calculate the mean values of (m s , m d , m v ) in each superpixel region to form representative points.

3.
Calculate similarity matrix Z between representative points and Laplacian matrix L rw .

4.
Calculate eigenvectors corresponding to the first K largest eigenvalues of matrix L rw as a matrix Q = [q 1 , q 2 , . . . q k ].
Cluster the row vectors of normalization matrix V with k-means.
The similarity matrix between representative points is calculated with LSC, which constructs a bipartite graph to decrease the memory and increase the speed of PolSAR image processing. For natural sandy land of complex characteristics, such as a large area, irregular objects [38,39], various features and intricate structures, LSC can be performed efficiently and flexibility. Moreover, the spatial information is effectively loaded with the help of superpixel segmentation algorithm, which further improves the final classification effect.

Experimental Results
As described in Section 2, there are various and complex types of features in the sandy land. Therefore, the area in the red color box is taken as the study area. The area is 1000 × 3000 pixels, mainly including residents (RT), roads (RD), semi-vegetation sand (SS), sandy land (DL), saline land (SL), vegetation (V), and lakes (L). The Pauli RGB image after Lee filtering is shown in Figure 5a. Test samples without overlap (i.e., from different fields) are randomly selected, as shown in Figure 5b, to evaluate the classification methods. The detailed information is shown in Table 2.

Decomposition Results
Three decomposition approaches are adopted to study area decomposition, and in-

Decomposition Results
Three decomposition approaches are adopted to study area decomposition, and investigated to determine which is more suitable. The first strategy is to adopt original hybrid Freeman/eigenvalue decomposition (HFED) [7]. The second is to utilize HFED with extended volume scattering (HFEDE), which is surely more reasonable than the first strategy. The extend volume scattering are applied to the ensemble average of dipole scatters or dihedral (horizontal or vertical) corner reflectors. The third approach is to use the ND, as described in Section 3. In 3.2.1 of Section 3, the cross-polarization scattering components of the PolSAR data are removed by OAC and PAR, which solves the problem that the arrangement direction of the building is not parallel to the radar azimuth direction. In 3.2.2 of Section 3, volume scattering mainly comes from the dihedral structure of the artificial target for more accurate modeling in artificial areas, while it comes from a cloud of randomly oriented dipole scatter such as vegetation in natural areas. The decomposition results obtained under the above three approaches are shown in Figure 6. At a glance, it is difficult to see the differences of the three images. Therefore, the three decomposition approaches corresponding to the percentage of the power values of the three polarization scattering powers (m s , m d , m v ) are shown in Figure 7, from which it can be seen that the volume scattering power ratio decreases from 52% with HFED to 41% using ND as proposed in this paper. At the same time, the ratio of surface scattering increases to 45%, which accords with the actual scattering characteristics of the region, and double-bounce scattering increases from 14% with HFED to 16% with HFEDE. Finally, the results of the decomposition approach more accurately reflect the scattering characteristics of real objects.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 23   For further analysis, three zones are selected and marked by red rectangles in Figure 5a and labeled Zone 1, Zone 2, and Zone 3. Each zone is 100 × 240 pixels, and the types of ground truth for these zones are lake areas, residents, and trees, respectively. It is well known that the dominant scattering powers of these three zones are surface, doublebounce, and volume, respectively, as shown in Table 3. In Zone 1, it can be seen that mean_m s as given by ND is 0.988, which is about 13.6% higher than for the other two approaches. In residential areas, because buildings and the ground constitute a large number of dihedral structures, double bounce scattering is dominant. In Zone 2, ND again performs best. Specifically, mean_m d from the proposed method is 0.604, which exceeds the mean values of HFED and HFEDE by 16.6% and 10.4%, respectively. However, in Zone 3, mean _ mv with HFED is the largest. In the ND, the cross-polarization power contributes to all of the three scattering models, which leads to the surface scattering power and double-bounce scattering power being higher than those generated by the current hybrid Freeman/eigenvalue decomposition techniques. On the whole, ND obtains better decomposition performance than the other methods.

Superpixels Generation Results
Two adjustable parameters related to the ASLIC algorithm are introduced in 3.3 of Section 3. The parameter β is used to balance the distance between the pixel value and the distance in the image plane space, and it usually has a fixed value of 1. N s is used to control the superpixel size of a block. This determines the similarity matrix Z, which affects the accuracy of the final classification result. At the same time, to produce better superpixel segmentation results, we test the segmentation effect when N s is 10, 30, and 50. The superpixel generation maps are shown in Figure 8. From the magnified view of the region in the box, we can see that ALSIC shows better boundary-preserving ability than SLIC for superpixel generation of PolSAR images, especially for some narrow or smaller targets. Specifically, the structures of target can obviously be observed after superpixels segmentation works with N s = 30, including good detail in the island and building areas. Compared with Figure 8b,e shows better segmentation quality, as the sizes of superpixels are similar, the borders of superpixels are smoother, and the edges of the island and building areas are preserved well. So, ASLIC outperforms SLIC for superpixel segmentation results under the same parameter N s .
To illustrate the superiority of ASLIC in superpixel generation of PolSAR image, we compare the performances of the ASLIC with generalized mean shift (MS) [30] and Ncut [30], and evaluate the performances of the superpixels by boundary recall (BR) [27] and achievable segmentation accuracy (ASA) [27]. BR represents the degree of maintains between superpixels and real boundaries. ASA indicates the degree of coincidence between superpixels and real targets. Figure 9 shows the ASA scores and the BR scores. As the number of superpixels increases, the performance of each algorithm is basically improved until it is relatively stable. Due to the complex terrain type, the gap between the ASA score and the BR score of these three methods are large to each other. However, ASLIC has the highest ASA score and BR score. It improves the BR score at 4000 superpixels to 95.0%, which is 17.6% and 18.7% higher than GMS and Ncut, respectively. At same time, the ASA score of ASLIC is the highest. Hence, we can conclude that the superpixel by ASLIC algorithm is suitable for the study area of Hunshandake Sandy Land.

Classification Results
In this section, the classification results corresponding to three methods are shown in Figure 10, which includes the results of unsupervised classification based on HFED and spectral clustering with superpixels of [33] (HED-SC), unsupervised classification based on the use of Neumann decomposition and Random Forest Classifier (ND-RF), and the proposed classification method (ND-LSC). We adopt the overall accuracy (OA) and Kappa coefficient to evaluate the performance of different methods. Accuracy is calculated by Equation (30):

Classification Results
In this section, the classification results corresponding to three methods are shown in Figure 10, which includes the results of unsupervised classification based on HFED and spectral clustering with superpixels of [33] (HED-SC), unsupervised classification based on the use of Neumann decomposition and Random Forest Classifier (ND-RF), and the proposed classification method (ND-LSC). We adopt the overall accuracy (OA) and Kappa coefficient to evaluate the performance of different methods. Accuracy is calculated by Equation (30): where TN is the number of correct classifications for all the training samples and NG is the total number of pixels in ground truth, P e represents the proportion of misinterpretation caused by accidental factors in the interpretation of remote sensing images. For each class, user's accuracies (UA) and producer's accuracies (PA) explained in Equation (31) are also calculated [27]. PA = TP TP+FP , U A = TP TP+FN (31) where TP is the true positive or the number of correctly detected class points, FP is the false positive or the number of incorrectly detected class points, and FN is the false negative. The classification accuracies of these three classification methods are listed in Table 4. It can be observed that the proposed method outperforms the others. On the one hand, the OA of the proposed method is 95.22%, while the values of the other two methods are 89.68% and 90.02%, respectively. ND-LSC has the best kappa coefficient (0.9404). On the other hand, for six out of seven classes, the proposed method performs better than ND-RF. For vegetation, the classification precision improves up to from 79.89% for HED-SC to 97.87% for the proposed method. In addition, for each class, the proposed method obtains better UA. In summary, our approach performs better on complex terrain classification. The findings in [11] also show that the parameters from the various polarimetric decompositions benefit the overall accuracy of classification methods. Polarimetric decompositions represent the scattering mechanisms of natural environments by scattering models. For the dominant scattering mechanism from different decomposition, the LSC scheme is able to better pick out small differences in the signal that characterizes the different classes. The violin plots introduced in Figure 11 illustrate three scattering mechanisms based on HFED, and ND decomposition shows differentiation between classes. Where Figure 11a,b show the distribution of the double-bounce scattering power coefficient m d found within each class, Figure 11c,d show the distribution of the surface scattering power coefficient m s found within each class Figure 11e,f show the distribution of the volume scattering power coefficient m v found within each class. It can be seen that the m d values from the HFED shows little differentiation between classes compared from the ND, which is generally associated with the double-bounce mechanism. However, the classes with the lowest UA and PA are the SL, which have significant overlap in double-bounce scattering, as seen in Figure 11a,b. In addition, for distribution of power coefficients m v and m s , compared with ND decomposition, the scattering mechanism based on HFED shows that the differences between classes are relatively small. In summary, the polarization parameters from ND benefit to the better performance of LSC classification in the sandy land.
To further illustrate the effectiveness of the proposed approach, six classification strategies are implemented to discuss the impact of accuracy indices (%) with different superpixels sizes. The parameter N s ranges from 10 to 50, with an interval of 10. Figure 12 displays the experiment results. It can be seen that almost every classification strategy has the highest accuracy when N s = 30. The classification accuracy of ND-SLIC-LSC exceeds those of HFED-SLIC-LSC and HFEDE-SLIC-LSC, because ND more accurately obtains the polarization scattering power (m s , m d , m v ). HFED-ASLIC-LSC has better classification accuracy than HFED-SLIC-LSC; because ASLIC generates superpixels, it has a stronger adaptation ability, and selection of spatial features is acquired to greatly improve classification accuracy. Thus, the proposed ND-ASLIC-LSC has the best accuracy-more than 94%.
Based on above discussions, several conclusions can be draw from experimental results. On the one hand, the proposed method (ND_LSC) obtains the best performance seen as OA (95.41%) and Kappa coefficient (0.943) in two categories of experiments, ie, it can better obtain a higher OA and Kappa coefficient by taking full advantage of the polarimetric features and spatial information. On the other hand, the performance using ND strategy is better than other decomposition strategies resulting from scattering powers (m s , m d and m v ) can be extracted more efficiently through the ND strategy. Secondly, the application of LSC is more in line with characteristics of its large area. The feature vector X combined with spatial coordinate information constitutes the input vector of ASLIC, which is used to generate superpixels.   To fully understand the consistency of classification results with real features, a field survey was conducted with Sangendalai town as the center. Several major types of land cover are selected for further analysis. Real photographs of land cover are shown in Figure 13, including vegetable I, vegetable II, semi-vegetation sand, lake I, lake II, saline land I, saline land II, sandy land I, sandy land II, road I, road II, and residents. The ground truth is used as a reference for the analysis of classification results. The geographic locations corresponding to the latitude and longitude of different land covers are recorded and displayed in Table 5. We compare the actual land cover with the result of the proposed classification method. It can be seen that Figure 13a,b show both types of vegetation, which correspond to two test samples in the green area. The green dense area represents vegetation II, and the sparse area represents vegetation I in Figure 10c. Most of the salinealkali land is distributed around the lake by observing Figure 13d-g. Correspondingly, the yellow area surrounds the purple in Figure 10c. There are two typical roads in Figure 13j,k, one surrounded by buildings and the other by vegetation. It can be seen from Figure 10c that some light blue thin lines cross areas and surround light green areas. Some light-blue blocks are seen by the lake. This phenomenon will be discussed below.

Discussion
We have proposed an unsupervised classification method (ND-LSC), which combines the advantages of the ND strategy, the ASLIC algorithm, and the LSC algorithm. The experimental results in Hunshandake Sandy land manifested that the proposed method accomplishes the best performances with regard to OA and Kappa. This is mainly attributed to two aspects. On the one hand, an ND strategy method is adopted to extract polarimetric scattering powers (m s , m d and m v ), which contribute to superpixels being more reasonably generated by the ASLIC algorithm. On the other hand, the unsupervised classificationbased LSC algorithm calculates the similarity matrix between representative points of superpixels generated by the ASLIC algorithm, which constructs a bipartite graph to solve the problem of insufficient memory and slow speed for PolSAR image processing. The following points in the experimental results are worth discussing.
• ND strategy can more effectively obtain polarization parameters than other decomposition strateges. There are two main reasons for this phenomenon. One is that the traditional decomposition strategy has a better performance for PolSAR imagery in other fields, but the characteristics of Hunshandake Sandy land are not very ideal. Another reason is that OAC and PAR are not performed and volume scattering models cannot adapt to the environment in natural and artificial areas during decomposition, which affects the decomposition result. After OAC and PAR strategy, which enhances the T 11 and T 22 elements of coherence matrix T 3 power. Among them element T 11 is relevant to surface scattering mechanism; element T 22 is relevant to double-scattering mechanism. In other words, classes dominated by surface scattering and double-scattering mechanism can be classify accurately. Therefore, the polarization parameters extracted by ND strategy can be utilized to generate superpixels in Hunshandake Sandy land.

•
There are the details that cannot be ignored of the classification result in Figure 10b, the obvious thing that can be observed is the road area around the lake. This phenomenon is not commission errors, classifying areas that are not roads as roads. This could be put down to the fact that both classes are predominantly surface scattering (water surface and roads surface) with relatively low surface roughness. The ground is very smooth after some small lakes degraded, which makes very similar to scattering mechanism of road. This would indicate that this specific environment has a polarimetric signature that could be associated with many classes.
However, there still are some disadvantages in the proposed method: (1) the algorithm structure is relatively complex and time-consuming when compared with conventional unsupervised classification methods. (2) This method has relative limitations and currently has a good classification result in the Hunshandake sandy land.

Conclusions
We have proposed an unsupervised classification method for Hunshandake Sandy Land. It integrates three key steps: the new decomposition (ND), superpixels segmentation, and the LSC algorithm. We have analyzed the capabilities of these steps to improve classification performance. Specifically, three classification methods and six classification strategies are discussed and applied to evaluate their effectiveness and feasibility. The analysis shows that the proposed method has superior classification performance on OA (95.41%) and kappa coefficient (0.943). Three years of data from Hunshandake Sandy Land are utilized to verify the superiority of the method. We conclude the following: (1) the ND strategy can better extract polarization scattering power compared with traditional decomposition, especially to distinguish sandy land types and (2) the ASLIC algorithm can further organize these ND-based polarization features combined with spatial coordinate information to generate superpixels, which can improve classification performance. Nevertheless, further investigation is still necessary for Hunshandake Sandy Land.