# Assessing the Long-Term Evolution of Abandoned Salinized Farmland via Temporal Remote Sensing Data
## 3 Methods
### Analyzing Temporal Remotely Sensed Imageries
#### 3.2.1 Land Cover Classification via _K_-Means Clustering  
_K_-means clustering [34], an unsupervised method that accepts arbitrary shapes of input for an individual year, was used to preliminarily classify Earth's surface. The applicability of unsupervised methods against supervised methods in different situations is discussed in Section 5 This algorithm assumes that the samples have \(k\) cluster centers and optimizes the cluster centers by minimizing the within-cluster distances (objective function) as follows:  
\[\text{obj}=\text{min}\sum\limits_{i=1}^{k}\sum\limits_{x\in C_{i}}||x-\mu_{i}|| ^{2} \tag{7}\]\[\mu_{i}=\frac{1}{C_{i}}\sum\limits_{x\in C_{i}}x \tag{8}\]  
where \(C_{i}\) is the cluster group and \(\mu_{i}\) is the cluster center for Group \(i\). The initial positions in the centroids were randomly selected from the samples. The variation in initial centroids introduced uncertainty into the final result. Therefore, the \(K\)-means algorithm was run 10 times with different centroid seeds in this study, and the best case of initial centroids was selected based on its rate of convergence (\(K\)-means clustering was deployed by the Python Sklearn package). After clustering, representative final clusters corresponding to main land types were manually selected to analyze their typical VI patterns (Figure 3) and to build a plant shape model to distinguish specific vegetation types (Section 3.2.2). The temporal reflectance dynamics of distinct surfaces are diverse due to the significant differences in the spectral reflectance characteristics of different landscapes. For a planting area, the number of peaks of the VI trajectory was determined by the number of cropping cycles. There is only one growing season in the Hetao Irrigation District due to the extremely low temperature in the winter. As a result, the VI dynamics of a cropping area show a single peak at the maximum leaf area index (LAI) stage (for instance, spring wheat reaches the maximum greenness at the heading stage at approximately DOY 160, while summer maize and sunflower peak at approximately DOY 210). The bodies of water show a negative temporal NDVI throughout the year that makes them easy to identify. The water body area accounts for a small proportion of the total area of the study area, and the spatial variability of water quality is relatively large. As a result, the standard error of NDVI is relatively larger for the body of water. The non-vegetative land has an NDVI with low variability, and the forest shows a long duration in growing season. However, the final classification results were unstable for different random seeds, although the initial centroids were optimized. For this reason, a large \(k\), which is far greater than the actual classes, was used to explore more potential land types. Then the threshold-based method was used to detect and merge the non-vegetative area and water body, and the \"shape-model-fitting\" method was used to detect and merge specific vegetation types from the final clusters. A higher \(k\) allows the algorithm to split the vegetative and non-vegetative areas into land-use types at a fine scale but at the cost of heavy computation. The value of \(k\) was determined by the trial-and-error method, and \(k\) = 20 was finally applied to make the uncertainty of the final detection of land type small. Then, the types of interest were identified and merged by predefined rules. Specifically, cluster centers with an annual average NDVI (average value of NDVI time series for each year) ranged from 0 to 0.2 were regarded as non-vegetative land, and water bodies were confirmed if the value was negative. For the vegetative area, a shape-model-fitting method is introduced in the next section to identify specific crops based on phenological metrics.  
Figure 3: Examples of cluster centers with annual temporal NDVI in the Hetao Irrigation District in 2019. The shaded width area denotes the standard error (calculated by the NDVI value of all pixels in this class) for each group.