*2.3. Methods*

The workflow used in this research is provided in Figure 2.

**Figure 2.** The workflow of crop type mapping by the integration of time series Sentinel-1 and Sentinel-2 features.

Firstly, a number of SAR features were extracted from time series Sentinel-1 data. It has been reported intensity and derived ratios (VV versus VH ratio, and the normalized ratio procedure between bands (NRPB)) can monitor the vegetation changes. Some studies indicated that SAR interferometric coherence can potentially improve the discrimination capability of different land cover types [25–28]. Thus, the backscatter intensity, interferometric coherence, and their derivative products were derived as SAR features for crop mapping. A by-product from interferometric pairs was also included; that is, the master versus slave intensity ratio computed respectively for the VH and VV polarization modes. In addition, we derived a statistic of time series SAR intensity, i.e., the amplitude dispersion computed respectively for VH and VV polarization. The amplitude dispersion was originally used in multi-temporal interferometric SAR (InSAR) for the initial selection of persistent scatters, which is a good indicator to distinguish vegetated surface and man-made structures, bare rocks, etc.

Secondly, optical features were derived from Sentinel-2 images. In addition to the conventional features, such as the visible, NIR (Near-infrared), SWIR (Short-wave infrared) bands, normalized difference vegetation index (NDVI), and normalized difference water index (NDWI), the red-edge spectral bands were also extracted. A total of 11 red-edge indices (will be detailed in Table 7) were calculated for each S2 acquisition, in order to thoroughly explore the contribution of red-edge wavelengths in crop type discrimination.

To deal with the high dimensional input features for the classifier, we proposed a recursive feature increment approach to select the optimal combination of SAR and optical features (detailed in Section 2.3.2) based on the random forest feature importance ranking.

A two-step hierarchical cotton mapping strategy (Figure 2) was implemented in this research. In step 1, the entire field site was classified into six land cover types, i.e., cropland, forest, desert, urban area, water body, and wetlands, so as to create a cropland mask. In the second step, the cropland was mapped into five different crop types, including chili pepper, corn, cotton, pear, and tomato. The automated feature selection method and RF classification were used both in step 1 and step 2, as illustrated in Figure 2.

#### 2.3.1. Data Processing and Feature Extraction

#### 1. Sentinel-1 Data Processing and Feature Extraction

All Sentinel-1 IW SLC images were pre-processed by the SNAP—ESA Sentinel Application Platform v7.0.0 (http://step.esa.int). Processing steps include co-registration, de-burst, and subset. To derived InSAR products, in total, 17 interferometric pairs were formed for each polarization mode using the acquisition on 26 March 2018 as the common master image, and all the subsequent acquisitions as the slave images. The next step is the de-speckling of intensity and accurate estimation of coherence. As pointed out by [29], using regular-shaped windows, conventional de-speckling methods average the values of neighboring pixels indiscriminately, leading to degraded image resolution and blurred edges between objects of different scattering characteristics. Furthermore, conventionally, the coherence is calculated by a regular-shaped sliding window, which unavoidably averages pixel values from different distributions (e.g., belonging to different land cover types), resulting in a biased estimation. In this study, a SHP DSI algorithm [23,24] was applied to the single look complex SAR data prior to SAR feature extraction. Here, we briefly summarize the principles and processing steps of the SHP DSI algorithm:

In SAR images, there are large numbers of distributed scatters (DSs) that exist in a resolution cell in cropland, forest, desert, etc. DSs cannot dominate the backscatter characteristics of a resolution cell, and the pixel appears as a random variable along the time dimension. When the SAR time series is sufficiently long, according to the central limit theorem, the vector sum of all the distributed scatters from a pixel can be defined as a complex Gaussian random variable. In a SAR image, we can identify if a pixel has the same statistical distribution with the other using the confidence interval. In this way, pixels can be divided into different statistically homogeneous pixel (SHP) subsets. The SHP DSI algorithm applies a modified Lee filtering method to the diagonal elements of the complex covariance matrix of an arbitrary pixel, on the basis of the SHP subset that the pixel belongs to, to obtain filtered time series intensity.

As for the interferometric coherence, for an arbitrary pixel *i*, the coherence coefficient is calculated as:

$$\gamma = \frac{\left| \sum\_{i=1}^{K} \mathbb{S}\_1(i)\mathbb{S}\_2^\*(i) \right|}{\sqrt{\sum\_{i=1}^{K} \left| \mathbb{S}\_1(i) \right|^2 \sum\_{i=1}^{K} \left| \mathbb{S}\_2(i) \right|^2}} \tag{1}$$

where *S*1(*i*) and *S*2(*i*) represent the corresponding complex values of the pixel *i* in the master and slave images, and *K* is the number of pixels in an SHP subset where the pixel *i* is located. Firstly, the coherence coefficient is calculated on the basis of each SHP subset for the purpose of removing the errors caused by the image texture. Secondly, a maximum-likelihood fringe rate algorithm [24] is applied to compensate for the interferometric fringe pattern. Finally, the bias in the coherence estimation is further mitigated using the second kind of statistical estimator proposed by [30].

After the implementation of the SHP DSI algorithm, we obtained the de-speckled backscatter intensity and bias-corrected coherence coefficient for each acquisition date. The units of *VH* and *VV* intensity were both converted decibels (dB). On this basis, a number of other Sentinel-1 features were generated as follows.

As one of the InSAR products, the master versus slave intensity ratio of the *VH* and *VV* polarization modes was calculated as

$$\begin{array}{l}\text{Intensity\\_VH\\_mst/slv\\_ratio} = \frac{\text{int}\_{V\bar{H},mt}}{\text{int}\_{V\bar{H},sh}}\\\text{Intensity\\_VV\\_mst/slv\\_ratio} = \frac{\text{int}\_{VV\bar{J},mt}}{\text{int}\_{VV\bar{J}v}}\end{array} \tag{2}$$

where *intVH*,*mst*, *intVV*,*mst* and *intVH*,*slv*, *intVV*,*slv* are respectively the backscatter intensity (in the unit of dB) of the master and slave images of each interferometric pair with the *VH*, *VV* polarization mode.

A statistic previously used in multi-temporal InSAR, i.e., amplitude dispersion is calculated as

$$amp\\_disp\_{VW} = \frac{std(amp\_{VW})}{mcan(amp\_{VW})}, \; amp\\_disp\_{VV} = \frac{std(amp\_{VV})}{mcan(amp\_{VV})} \tag{3}$$

where *ampVH* and *ampVV* represent the amplitude of each Sentinel-1 acquisition of the *VH* and *VV* polarization mode.

A by-product of backscatter intensity, i.e., *VV* versus *VH* intensity ratio is calculated as

$$Intensity\\_VV/VH\\_ratio = \frac{int\_{VV}}{int\_{VH}}\tag{4}$$

where *intVV* and *intVH* are respectively the backscatter intensity (in the unit of dB) of each Sentinel-1 acquisition.

Another intensity-based ratio, the normalized ratio procedure between bands (NRPB) [31] is calculated as

$$NRPB = \frac{\textit{int}\_{VH} - \textit{int}\_{VV}}{\textit{int}\_{VH} + \textit{int}\_{VV}} \tag{5}$$


**Table 7.** Spectral indices calculated from Sentinel-2 data.

In summary, we derived 10 subcategories of features from Sentinel-1 data. Adding features of different acquisition dates together, 142 Sentinel-1 features were produced, geocoded, and resampled to 10 m resolution.

#### 2. Sentinel-2 Data Processing and Feature Extraction

Sentinel-2 data were firstly processed from Level-1C to Level-2A surface reflectance using a Sen2Cor—ESA Sentinel-2 Level 2A processor (https://step.esa.int/main/third-party-plugins-2/sen2cor/ (accessed on 30 October 2019)). Cirrus cloud correction was done during the Level-1C to Level-2A processing by Sen2Cor. Pixels covered by thick clouds were masked out using the opaque clouds band in the Level-2A product metadata. All spectral bands of 10 m and 20 m spatial resolution, including the red-edge bands (Band 5, Band 6, and Band 7) were included in the analysis. The 20 m resolution bands were resampled to 10 m resolution for further processing. Apart from the normalized difference vegetation index (NDVI) and normalized difference water index (NDWI), 11 red-edge indices (Table 7) were generated for each Sentinel-2 acquisition.

All geocoded Sentinel-1 features were aligned with the Sentinel-2 features. In summary, 33 subcategories, a total of 326 Sentinel-1 and Sentinel-2 features, were produced for subsequent analysis.

#### 2.3.2. Optimal Feature Selection and Classification

#### 1. Random Forest Classification

The random forest (RF) classifier [36] was deployed in this study for both cropland extraction and crop type classification. Random forest is a machine learning algorithm for classification, which consists of an ensemble of decision trees where each tree has a unit vote for the most popular class to classify. RF has been widely used for remote sensing classification, as it runs efficiently on large databases and performs well with high-dimensional features.

In this study, we used the RF classifier from the Scikit-learn Python module [37] for classification. The two key parameters, i.e., the number of trees and the number of features in each node, were chosen by analyzing the out-of-bag (OOB) errors. As a result, the number of trees was set as 450 in step 1 and 550 in step 2, and the number of maximum features in each tree was set as the square root of the total number of input features in both steps.

#### 2. Optimal S1 and S2 Input Variable Selection Using Recursive Feature Increment (RFI) Method

Inspired by the backward feature elimination method [22], we proposed a forward feature selection approach to obtain the optimal combination of S1 and S2 features, which is referred to as the recursive feature increment (RFI) method.

Firstly, a feature importance ranking was obtained by the permutation importance of random forest calculated for each feature. The permutation importance assesses the impact of each feature on the accuracy of the random forest model. The general idea is to randomly permute the feature values and measure the decrease of the accuracy due to the permutation. In this study, we utilized the permutation importance function from the ELI5 package to estimate the feature importance.

Secondly, using the feature ranking with the most important feature placed at the top, RF classification was implemented by recursively considering bigger and bigger feature sets, starting from one feature and adding one at a time. For the feature set at each iteration, an RF classifier is parameterized and assessed using stratified k-fold cross-validation, with the mean overall accuracy (OA), mean kappa coefficient, and mean f1 score of all classes, as well as the f1 score of individual class recorded.

Finally, by analyzing the time series of the accuracy metrics, for example, the kappa coefficient, the feature set corresponding to the iteration achieving the highest accuracy is selected as the optimal combination of S1 and S2 features.

#### 3. Accuracy Assessment

The performance of classification was assessed using stratified k-fold cross-validation. The entire sample dataset was split into k stratified "folds". The folds were made by preserving the percentage of samples for each class. In the case of unbalanced samples, this is to ensure that each fold is a good representative of the whole. For each fold in the sample dataset, the classification model was trained on k-1 folds and tested on the kth fold. This was repeated until each fold had served as the test set. Three accuracy metrics were used in this study, i.e., the overall accuracy (OA), Cohen's kappa coefficient, and F1 score of each class. All these metrics were averaged over the k folds.

The OA is calculated from the confusion matrix by adding the number of correctly classified sites together and dividing it by the total number of the reference sites. It is to test what proportion were mapped correctly out of all of the reference sites. From the confusion matrix, the producer's accuracy (PA, also referred to as recall) and the user's accuracy (UA, also referred to as precision) can also be calculated. UA (precision) is the ratio of correctly predicted positive observations to the total predicted positive observations; PA (recall) is the ratio of correctly predicted positive observations to all the observations in the actual class. The kappa coefficient is to evaluate if the classification does better than randomly assigning values by taking into account the possibility of the agreement occurring by chance. In the case of uneven class distribution, the F1 score is a more useful metric than OA, as it is a weighted average of precision and recall. The F1 score for each class is calculated as follows:

$$F\_1 = 2\frac{precision \times recall}{precision + recall} \tag{6}$$

#### **3. Results**

#### *3.1. Step 1: Cropland Extraction*

A random forest (RF) classifier was applied to the Sentinel-1 and Sentinel-2 features for land cover classification, so as to generate a cropland mask. As the purpose of this step is to extract a cropland mask, we use the F1 score of cropland (Figure 3) as the accuracy metric for the RFI feature selection (Section 2.3.2). The top 114 features in the importance ranking were selected. The feature importance scores of the top six features are displayed in Figure 4, with corresponding boxplots of different land cover types shown in Figure 5.

**Figure 3.** The F1 score of cropland recorded in each iteration of the recursive feature increment (RFI) process, reaching the maximum at the 114th iteration. Therefore, the top 114 features in the importance ranking were chosen as the optimal feature set for land cover classification in step 1.

**Figure 4.** The feature importance scores of the top six features selected by the RFI method for cropland extraction, in descending order. Each importance score was normalized and converted to percentage. Each feature is named by its subcategory and acquisition time in the format of "yyyymmdd".

**Figure 5.** Boxplots of the top six most important features of different land cover types. Each feature is named by its subcategory and acquisition time in the format of "yyyymmdd".

In Figure 5, the top six features exhibit a good ability to distinguish different cover types, which reflects the reliability of the RFI feature selection method.

The random forest classifier was parameterized using the selected 114 features to generate a land cover map (Figure 6a). Then, a cropland mask (Figure 6b) was produced from the land cover classification results. The land cover classification accuracy was assessed using 10-fold cross-validation (Table 8).

**Figure 6.** (**a**) Land cover classification map; (**b**) cropland mask extracted from the land cover classification results.



#### *3.2. Step 2: Crop Type Mapping*

#### 3.2.1. Optimal S1 and S2 Feature Combination and Crop Type Mapping

Using the RFI feature selection method, we obtained the optimal combination of Sentinel-1 and 2 features for crop type discrimination. In this case, the kappa coefficient was used as the main accuracy metric to determine the final feature set. The mean OA achieved its maximum at the same time as the kappa coefficient, at the 113th iteration. Thus, in total, 113 features were selected. The mean OA and mean kappa coefficient averaged over the k-fold cross-validation are plotted in Figure 7. The feature importance scores of the top six features selected for crop type classification are shown in Figure 8.

**Figure 7.** The mean overall accuracy (OA) and mean kappa coefficient of fivefold cross-validation recorded for each iteration during the RFI feature selection process.

**Figure 8.** The feature importance score of the top six features selected by the RFI method for crop type classification, in descending order. Each importance score was normalized and converted to a percentage. Each feature is named by its subcategory and acquisition time in the format of "yyyymmdd".

The boxplots of the top six features (Figure 9) reveal the good separability of different crop types, which proves the rationality of the RFI feature selection results. For example, the intensity of the VH mode reflects a good ability to distinguish pear from other crop types in March. Pear values in Band 8 (NIR) in July show a clear distinction from other crops. Corn can be easily separated in the red-edge index NDVIre2n in September. Cotton is clearly distinguished in Band 11 (SWIR 1) in June.

Based on the optimal combination of S1 and S2 features, a random forest classification model was trained using the ground samples and applied to the whole study area to produce a crop distribution map (Figure 10). The classification accuracy was assessed by fivefold cross-validation (Table 9). Besides, the accuracy was also verified by a "one sample per field" method; that is, from an individual validation field (as listed in Table 3), only one sample was extracted. In this method, the validation samples are independent from each other. The corresponding accuracy metrics are listed in Table 10.

**Figure 9.** Boxplots of different crop types of the top six important features selected by the RFI approach. Each feature is named by its subcategory and acquisition time in the format of "yyyymmdd". (**a**) Boxplots of the band 11 values from the 3 June 2018 acquisition of Sentinel-2; (**b**) boxplots of the band 6 values from the 17 August 2018 acquisition of Sentinel-2; (**c**) boxplots of the NDVIre2n indices derived from the 1 September 2018 acquisition of Sentinel-2; (**d**) boxplots of the VH intensity values from the 26 March 2018 acquisition of Sentinel-1; (**e**) boxplots of the band 8 values from the 23 July 2018 acquisition of Sentinel-2; (**f**) boxplots of the band 7 values from the 17 August 2018 acquisition of Sentinel-2.

**Figure 10.** Crop distribution map of the study area derived from random forest (RF) classification using the combined Sentinel-1 and Sentinel-2 feature set.


**Table 9.** Accuracy of step 2 crop type classification using the optimal combination of S1 and S2 features, assessed by stratified fivefold cross-validation using the ground samples.

**Table 10.** Accuracy of step 2 crop type classification using the optimal combination of S1 and S2 features, assessed by stratified fivefold cross-validation using one sample from each validation field.


From Figure 10, we can see that pears are mainly planted in Korla County; the four counties near Bosten Lake (Bohu, Yanqi, Hejing, and Hoxud) are the main production areas of chili peppers and tomatoes; cotton is mostly cultivated in Yuli County on the edge of the Tarim Basin. Corn cultivation spreads across the whole study area.

3.2.2. Performance Comparison of SAR Features Filtered by Different Methods in Crop Classification

In this section, we compare the crop classification results obtained using SAR features processed by conventional de-speckling methods and the SHP DSI method. The input variables include all the SAR features mentioned in Section 2.3.1. In the conventional procedure, the intensity and intensity-derived features were de-speckled using a 7 × 7 refined Lee filter, and the coherence was estimated using a 7 × 7 sliding window. In the SHP DSI method, statistically homogeneous pixels (SHPs) were identified based on a 5 × 9 window using the intensity stack formed by 18 acquisitions. Then, the de-speckled intensity and accurately estimated coherence were obtained by the procedure described in Section 2.3.1. An example intensity of a subset area extracted from the original SAR intensity, intensity filtered by refined Lee method, and intensity filtered by the SHP DSI method are compared in Figure 11. Figure 12 shows the example coherence of the same subset area extracted from the interferometric coherence estimated respectively using a 7 × 7 sliding window and the SHP DSI method. The crop classification accuracy obtained by using SAR features processed by different filtering algorithms is listed in Table 11.

**Figure 11.** A subset of VH intensity on 26 March 2018 extracted from (**a**) Original synthetic aperture radar (SAR) intensity; (**b**) SAR intensity filtered by the refined Lee method; (**c**) SAR intensity filtered by the SHP distributed scatterer interferometry (DSI) method.

**Figure 12.** A subset of VV coherence on 7 April 2019 extracted from (**a**) Coherence coefficient estimated by a 7 × 7 sliding window; (**b**) Coherence coefficient estimated by the SHP DSI method with bias mitigation.

**Table 11.** Accuracy metrics of crop type classification using SAR features processed by different filters. SHP: statistically homogeneous pixel.


In Figure 11, it is evident that the speckles in original SAR intensity are both reduced by a refined Lee filter and SHP DSI filter, but the sharpness and details of each image is better preserved in (c). In Figure 12, it is found the coherence over a water body is overestimated in (a). This bias is reduced in (b) by using the SHP DSI method. It is clear the coherence in (b) exhibits less salt and pepper-like noise.

In Table 11, it is demonstrated that the SAR features filtered by the SHP DSI method yield the best accuracy in crop type classification. Compared with the results of using a refined Lee algorithm, the overall accuracy is increased by 6.25%, the kappa coefficient is improved by 0.08, and the F1-scores of each crop type are all improved.

3.2.3. Performance Comparison of Features from Different Sources in Crop Type Mapping

A comparison was conducted on the performance of crop type classification using four groups of features. The first group includes only Sentinel-1 features; the second group contains only conventional optical features exclusive of red-edge features; the third group comprises all Sentinel-2 features inclusive of red-edge contribution; the last group consists of the both SAR and optical features. The same feature importance ranking obtained in Section 3.2.1 was used for all groups. For each feature group, the optimal feature set was individually determined by the RFI feature selection method, using the kappa coefficient of each group as the accuracy metric. The mean OA and kappa coefficient reached the maximum at the same iteration index for each test, as summarized in Table 12. The corresponding F1 scores of different crop types are compared in Figure 13.


**Table 12.** Accuracy assessment of crop type discrimination using different groups of features. The mean overall accuracy (OA) and kappa coefficient were averaged over fivefold cross-validation in both multi-sample per field validation and one sample per field validation.

**Figure 13.** Comparison of crop type mapping results using different feature combinations. Four groups of features were tested, the first group contains only SAR (Sentinel-1) features, the second group contains only optical (Sentinel-2) features without the red-edge contribution, the third group has all of the Sentinel-2 features, the fourth group includes both SAR and optical (Sentinel-1 and 2) features.

In both Table 12 and Figure 13, it is evident that the combination of SAR and optical features achieved the best accuracy in crop type discrimination. The inclusion of red-edge bands and indices improved the mean OA and kappa coefficient respectively by 1.84% and 0.03 (3.06% and 0.04% in multi-samples per field validation), compared with the results using other Sentinel-2 features. By integrating the optical and SAR features, the mean OA and kappa coefficient were improved by 1.58% (1.55% in multi-samples per field validation) and 0.02% compared with the results only using Sentinel-2 features. In Figure 13, among the five crop types, the F1 score of chili pepper and corn were significantly improved by the inclusion of red-edge features.

#### **4. Discussion**

#### *4.1. Importance of Features from Di*ff*erent Sources*

The accumulated importance scores were calculated by subcategories (Figure 14), different data sources (Figures 15a and 16a), and the month of acquisition (Figures 15b and 16b).

**Figure 14.** Accumulated feature importance scores of features selected for (**a**) step 1: cropland extraction and (**b**) step 2: crop type discrimination, calculated for each subcategory regardless of acquisition time. Each importance score was normalized and converted to a percentage. Each feature is named by its subcategory and acquisition time in the format of "yyyymmdd".

**Figure 15.** The accumulated importance scores of features selected in step 1: cropland extraction in three groups. The first group contains only Sentinel-1 features; the second group comprises only Sentinel-2 features exclusive of the red-edge features; the third group includes only the red-edge features. (**a**) Feature scores in the three groups calculated regardless of the acquisition time; (**b**) Feature scores in the three groups calculated for each month. The VH amplitude dispersion, as a single-phase feature, is plotted on the rightmost bar. Each importance score was normalized and converted to a percentage.

**Figure 16.** The accumulated importance scores of features selected in step 2: crop type discrimination in three groups. The first group contains only Sentinel-1 features; the second group comprises only Sentinel-2 features exclusive of the red-edge features; the third group includes only the red-edge features. (**a**) Feature scores in the three groups calculated regardless of the acquisition time; (**b**) Feature scores in the three groups calculated for each month. The VH amplitude dispersion, as a single-phase feature, is plotted on the rightmost bar. Each importance score was normalized and converted to a percentage.

Among the features selected for cropland extraction (Figure 14a), it is found that the SWIR 2 band (Band 12) contributed the most in the cropland extraction, followed by the Sentinel-1 VH intensity and Band 11 (SWIR 1), again the red-edge 1 (Band 5), and NDWI.

In the features selected for crop type discrimination (Figure 14b), the VH intensity exhibited the highest accumulated score for crop type classification, followed by the NDVI, Band 8 (NIR), Band 6 (red-edge 2), Band 4 (red), and Band 5 (red-edge 1). The optimal feature set comprises six subcategories of SAR features, 13 subcategories of red-edge features, and nine subcategories of conventional optical features. The contribution of features from different data sources will be explored in detail later.

It should be noted that the VH intensity held a significant share in both land cover classification and crop discrimination. The intensity ratio between bands, i.e., S1 NRPB and VV versus VH intensity ratio were only used in land cover classification. As for InSAR products, it is found that the VV coherence was selected in step 2. No coherence was selected in step 1. The coherence feature with the highest score in its subcategory, i.e., the VH coherence on 17 August 2018, ranked 118 in the RF permutation importance results in step 1. Its importance score was 0.27%, which is comparable to the score of 0.29% of the last chosen feature, NDre2 (red-edge) index on 23 July 2018. The amplitude dispersion of VH polarization, the master to slave intensity ratio in both polarization modes, were used in both step 1 and step 2.

All conventional optical features and the three red-edge spectral bands contributed to both the land cover classification. Regarding the red-edge indices, "NDVIre3" was only used in step 1; all the other 10 indices were selected for both cropland extraction and crop type identification.

As shown in Figures 15a and 16a, the three groups of features held similar proportions in both steps. The conventional optical features had the largest proportion in step 1 and step 2; the SAR features accounted for −15% in step 1 and −16% in step 2; the red-edge features accounted for −24% in step 1 and −22% in step 2.

Comparing Figures 15b and 16b, the features selected in step 1 and step 2 had significantly different temporal distribution. For land cover classification, the features in May had the largest proportion, followed by September, October, and June. Red-edge features contributed more in September, June, and July. In crop type classification, it is clear that features in August held the biggest share, of which red-edge features had a significant proportion. Compared with other months, the red-edge wavelengths performed the best in August, when most crops were ripening. This is likely associated with the sensitivity of the red-edge wavelengths to differences in the leaf structure and chlorophyll content of crops. Conventional optical features (visible, NIR, and SWIR wavelengths) were dominant from May to July and in October, and held a significant share in August and September.

In both steps, there were no optical features available in March and April, but some SAR features show good separability in the early season. SAR features in June and July were not selected for either of the steps. The proportion of SAR features re-increased in August and at the end of harvesting season in October. The amplitude dispersion was used in both steps, as a single-phase feature, plotted on the rightmost bar in both Figures 15b and 16b.

In addition, we examine the correlation between selected features in both step 1 and step 2, as shown in Figure 17. In both step 1 and step 2, most of the selected features showed a low correlation, as revealed by the histograms in (c) and (d).

#### *4.2. The Contribution of SAR Features in Crop Type Discrimination*

The electromagnetic radiation of the C-band cannot penetrate through vegetation cover. Therefore, the radiation from the Sentinel-1 C-band SAR sensor reflects the interaction between the radar signal and the ground surface cover, which explains the sensitivity of Sentinel-1 SAR to the vegetation cover.

**Figure 17.** Heat maps showing correlation between selected features. (**a**) Correlation heat map of selected features in step 1; (**b**) correlation heat map of selected features in step 2. In both (**a**,**b**), the feature indexes follow the feature rankings obtained through the RF feature importance score. Correlation close to 1 or −1 indicates high positive or negative correlation between features. The diagonal elements in both correlation matrices are self-correlation coefficients, so they constantly equal one. (**c**) The histogram of correlation between selected features in step 1; (**d**) the histogram of correlation between selected features in step 2. In both (**c**,**d**), the histograms were calculated after removing the diagonal elements of the correlation matrices and converting a negative correlation coefficient to corresponding positive values. Thus, the '0' in the histogram indicates low correlation, while '1' indicates high correlation.

It has been reported that VH polarization is more sensitive to the structural and geometrical arrangements of plants [11,38]. This is in accordance with the result that VH intensity bands held the biggest proportion in the selected SAR features (also in all features) for crop type discrimination (Figure 14b). We examined the selected SAR features and found that 13 out of 22 features were in VH polarization mode, including six features of the VH intensity (Figure 18) and six features of the master versus slave VH intensity ratio (Figure 19), as well as the VH amplitude dispersion (Figure 20a). Pear is the most distinguishable crop in most of these features. This is highly likely because SAR intensity is sensitive to target structure, as the planting density and vegetation height in pear fields are significantly different from those in other crop fields. It is interesting to notice that the VV coherence on 7 April 2018 revealed a good ability to differentiate cotton from other crop types (Figure 20b). In Figure 20b, we can see the VV coherence is significantly low in the cotton field. According to local crop phenological calendars (Section 2.1), cotton is the crop that is most likely to be sown in early April. Other annual crops are usually sown or transplanted in mid-April or later. Thus, this is possibly related to the changes of the cotton field surface in the sowing season (early April) from bare soil to plastic-mulched land, which leads to a fast decrease of coherence in 12 days, while the other crop fields almost remain as before. In addition, chili pepper can be easier recognized from the VH intensity on 16 October 2018 (Figure 18f).

**Figure 18.** Boxplots of selected VH intensity features of different crop types, spanning the time from 26 March 2018 to 13 May 2018, as well as a feature on 16 October 2018. Each feature is named by its subcategory and acquisition time in the format of "yyyymmdd". In (**a**–**e**), pear values show an apparent distinction from other crops, while in (**f**), chili can be easier recognized.

**Figure 19.** Boxplots of the top three scoring master versus slave VH intensity ratio of different crop types. Each feature is named by its subcategory and acquisition time in the format of "yyyymmdd". (**a**) Boxplots of the master vs. slave VH intensity ratio derived from the 29 August 2018 acquisition of Sentinel-1; (**b**) boxplots of the master vs. slave VH intensity ratio derived from 10 September 2018 acquisition of Sentinel-1; (**c**) boxplots of the master vs. slave VH intensity ratio derived from 4 October 2018 acquisition of Sentinel-1.

**Figure 20.** Boxplots of selected SAR features of different crop types. (**a**) Coherence coefficient in VV polarization mode; (**b**) amplitude dispersion in VH polarization mode. Each feature is named by its subcategory and acquisition time in the format of "yyyymmdd".
