Automatic Crop Classification Based on Optimized Spectral and Textural Indexes Considering Spatial Heterogeneity

Abstract

Crop recognition with high accuracy at a large scale is hampered by the spatial heterogeneity of crop growth characteristics under the complex influence of environmental conditions. With the aim to automatically realize large-scale crop classification with high accuracy, this study proposes an automatic crop classification strategy considering spatial heterogeneity (ACCSH) by combining the geographic detector technique, random forest average accuracy model, and random forest classification model. In ACCSH, spectral and textural indexes that can quantify crop growth characteristics and environmental variables with potential driving effects are first calculated. Next, an adaptive spatial heterogeneity mining method based on the geographic detector technique is proposed to mine spatial homogeneous zones adaptively with significant differentiation of crop growth characteristics. Subsequently, in view of the differences in crop growth characteristics and key classification indexes between spatial homogeneous zones, correlation analysis, and random forest average accuracy are combined to optimize classification indexes independently within each zone. Finally, random forest is used to classify the target crop in each spatial homogeneous zone separately. The proposed ACCSH is applied to automatically recognize crop types, specifically wheat and corn, in northern France. Results show that kappa coefficients of wheat and corn using ACCSH are 15% and 26% higher than those of classifications at the global scale, respectively. In addition, the index optimization strategy in ACCSH shows apparent superiority. Kappa coefficients of wheat and corn are 5–18% and 9–42% higher than those of classifications based on non-optimized indexes, respectively. In general, ACCSH can automatically realize crop classification with a high precision that suggests its reliability.

1. Introduction

Spatial distribution of crop planting is crucial to achieve accurate agricultural management and production guidance. Crop classification with high precision is in demand. The current widely used method is crop classification based on remote sensing image data, given its capacity to obtain the spatial distribution of a crop planting structure [1].

In most relevant studies, crop classification based on remote sensing data is carried out from a global perspective of the study area [2,3]. However, crop growth characteristics are highly related to increasing environmental conditions. Under the driving effects of multiple key environmental factors, crop growth characteristics inevitably exhibit strong spatial heterogeneity, especially at large scale, which can give rise to barriers for precise crop classification. Thus, considering spatial heterogeneity is an important research direction to improve crop classification accuracy.

Several researchers proposed layered strategies [4,5,6,7,8] to address spatial heterogeneity’s influence on crop classification, which can be generally divided into two types: methods based on climatic conditions and those based on management zones. Hao et al. [4] divided China into eight zones based on its geographic position and climate and established classification guidelines for each zone. The zoning method effectively increases crop classification accuracy. Aparecido et al. [6] divided their study area into sub-regions according to segmentation temperature and precipitation. NDVI index segments were used by Chen et al. [7] and Santos et al. [8] to create management zones for multilayer crop classification. The zoning classification strategy [4,5,6,7,8] can also successfully reduce the effect of spatial heterogeneity on crop classification. The accuracy of crop recognition based on spatially homogeneous zones is higher than that of the global scale. However, existing studies on the layered crop classification of zones generally use artificial division with sufficient prior knowledge, which are effective in certain situations, especially when researchers possess sufficient prior knowledge of the study area. However, at times, prior knowledge is difficult to obtain. Additionally, crop growth is highly relevant with complex multiple environmental variables, under which current approaches face difficulties to adaptively mine the spatial heterogeneity patterns and mechanism of crop growth characteristics. Thus, the effects of multiple environmental influences on crop growth necessitates an adaptive layered crop classification approach that considers geographical heterogeneity.

The geographic detector (GD) model proposed by Wang et al. [9] can effectively detect the key driving factors that cause certain spatial heterogeneities. The GD has been widely and successfully applied in many fields [9,10]. Meng et al. [11] used the GD to explore the internal driving factors of vegetation change in Inner Mongolia, effectively detecting and analyzing the driving effects of terrain, climate, annual precipitation, and other factors. Han et al. [12] successfully used a GD to explore the influence of environmental factors on soil-total phosphorus content in soil at different erosion levels. Given its established advantages and wide applicability in mining key factors and spatial homogeneous zones, the GD shows potential in detecting spatial heterogeneity characteristics and mechanisms of crop growth characteristics under the influence of multiple environmental variables.

Proper indexes are crucial for effective crop classification. Spectral features have proven effective and are widely used to measure crop growth characteristics [13,14,15]. In addition, several scholars carried out crop classification based on textural features that can effectively reflect crop characteristics and distinguish between crops [16,17,18]. The present study intends to integrate spectral and textural features as potential classification indexes. However, specific spectral or textural features in specific months reflect differences in key growth characteristics between crops and may vary with crop growth spatial homogeneous zones. Therefore, classification indexes must be optimized to further reduce redundancy and successfully improve accuracy. Index optimization strategies based on correlation and machine learning, such as the average accuracy reduction (RFAA) method, are confirmed to be effective in obtaining key indexes and improving crop identification accuracy. These strategies are used to optimize the indexes for crop classification in this study [19,20].

In the study area, we propose a novel adaptive crop clustering strategy (ACCSH) comprising three phases to achieve satisfactory crop classification considering the spatial heterogeneity of crop growth characteristics. In phase 1, an adaptive spatial heterogeneity mining method (SHGD) is proposed to mine spatial homogeneous zones of crop growth characteristics by combining a GD and clustering. In phase 2, classification indexes are separately optimized in each spatial homogeneous zone. On this basis, the classification model using random forest is constructed in phase 3. The remainder of this paper is organized as follows. Section 2 describes the study area and dataset. Section 3 elaborately illustrates the proposed algorithm and the corresponding accuracy analysis. Section 4 presents the experiments on real-world applications that are carried out to test the proposed method. Finally, the discussion and conclusion are presented in Section 5 and Section 6, respectively.

2. Study Area and Dataset

2.1. Study Area

The study area is the northern part of the French mainland, excluding Corsica (Figure 1), which is relatively flat. This area experiences warm weather and some rains all year round and is thus suitable for the growth of various crops. Agriculture is well developed with a high degree of modernization. The main crops widely grown are wheat, corn, and sugar beet. Since France is an important supplier and exporter of agricultural and sideline products globally, efficient and accurate monitoring of crops in this region is essential.

2.2. Data Acquisition and Analysis

The sample dataset with crop types and geographical location information is derived from the 2009 European LUCAS dataset (Figure 1). Wheat and corn are the two main crops concentrated in northern France and are therefore used as the study crop types. For the study area, 325 wheat samples, 600 non-wheat samples, 213 corn samples, and 545 non-corn samples were obtained. The remote sensing data are obtained from the 2009 MODIS13Q1 product (Figure 2). The Savitzky–Golay smoothing filter is used to reduce noise and retain effective spectral information. The environmental data (topography, temperature, precipitation) are obtained from WorldClim [21] (https://www.worldclim.org/data/index.html (accessed on 29 November 2009)), where the monthly average image data are selected for climate and precipitation. The spatial resolution is 2.5 min (~12 km²).

Phenology refers to the cyclical changes of crops in response to external environmental conditions. The varying phenological information and growth characteristics among crops are the key to their differentiation. The relevant phenological information on crops in the study area is obtained through a series of previous studies [22,23,24]. Figure 3 shows that the growth cycle of wheat is from October to July, and corn is from April to November. These months of the growth cycle must be included in the calculation of classification indexes.

3. Methodology

Spectral and textural indexes that can characterize crop features are discussed in Section 3.1. Environmental variables related to spatial heterogeneity are given in Section 3.2. The proposed ACCSH to classify crops is explained in Section 3.3. The classification evaluation indicators for assessing the classification effect are described in Section 3.4. The procedure is given in Figure 4.

3.1. Construction of Classification Indexes

3.1.1. Construction of Spectral Indexes

The MODIS13Q1 product captures images across five bands, including NDVI, red, blue, near-infrared (NIR), and middle-infrared (MIR), which were adopted as spectral indexes for 12 months in 2009. Additionally, different growth characteristics between crops are mainly reflected in the NDVI time series curves (Figure 3) in terms of spectral information [25,26,27]. Two quantitative indexes [26,27], curve trend $T_i$ and curve difference $D_i$, are proven to be effective to distinguish different crops. $T_i$ and $D_i$ are then calculated based on the NDVI values of crop phenology according to the equations in

Table 1. Spectral indexes and interpretation.

Indicator	Calculation Formula	Meanings
$n$-NDVI	-	The NDVI value in month $n$.
$n$-RED	-	The value of red band in month $n$.
$n$-BLUE	-	The value of blue band in month $n$.
n-NIR	-	The value of NIR band in month $n$.
n-MIR	-	The value of MIR band in month $n$.
$T_i\left(j\right)$	$T_i\left(j\right)=\frac{\sum (P_{in}-\overline{P_i})(\overline{T_{jn}}-\overline{T_j})}{\sqrt{\sum {(P_{in}-\overline{P_i})}^2\sum {(\overline{T_{jn}}-\overline{T_j})}^2}}$	Similarity between sample $i$ and crop type $j$ in NDVI index, $j$ is the target crop.
$D_i\left(j\right)$	$D_i\left(j\right)=\sqrt{\sum {(P_{in}-\overline{T_{jn}})}^2}$	Difference between sample $i$ and crop type $j$ in the NDVI index, $j$ is the target crop.

$i$ is the label of the classified sample; $\overline{T_j}$ is the average NDVI value of samples of crop type $j$ during the phenological period of crop type $j$. $T_{jn}$ represents the average NDVI value of the crop type $j$ in the month $n$. $n$ is the label of each month during the phenological period of crop type $j$. $\overline{P_i}$ is the average NDVI value of sample $i$ during the phenological period of crop type $j$, and $P_{in}$ represents the NDVI value of sample $i$ in the month $n$. Time phase number $n$ in this study is calculated using the phenological period of crop $j$.

. $T_i$ reflects the fluctuating trend of crop growth change and measures the correlation with the pseudo phenology curves [26]. High $T_i$ means more similarities between the sample $i$ and the target crop. $D_i$ measures the distance difference with the pseudo phenology curves and judges the difference between curves [27]. Low $D_i$ indicates a small difference between sample $i$ and the target crop. In summary, spectral bands, $T_i$ and $D_i$ were used as spectral indexes in this study.

3.1.2. Construction of Textural Indexes

As a commonly used texture algorithm, a gray-level co-occurrence matrix [28] was adopted to calculate textural indexes based on the images in Figure 2. Textural indexes include the mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, and correlation (

Table 2. Textural indexes and interpretation.

Index	Meaning
$i-B1$	The value of textural index mean in month $i$
$i-B2$	The value of textural index variance in month $i$
$i-B3$	The value of textural index homogeneity in month $i$
$i-B4$	The value of textural index homogeneity in month $i$
$i-B5$	The value of textural index dissimilarity in month $i$
$i-B6$	The value of textural index entropy in month $i$
$i-B7$	The value of textural index second moment in month $i$
$i-B8$	The value of textural index correlation in month $i$

$i$ is calculated using the phenological period of the target crop. When the target crop is wheat, $i=\{\mathrm{10,11,12,1}\cdots ,7$}; when the target crop is corn, $i=\{\mathrm{4,5},\cdots ,11\}$.

). The mean reflects the relationship between the brightness of the target pixel and its neighbors. Variance indicates the degree of grayscale variation in the local area of the image. Homogeneity shows the concentration of pixels in the matrix. The contrast reflects the degree of gray difference in the local area. Dissimilarity is the linear correlation of the contrast of pixels in the local area. Entropy was used to measure the textural information, and the richer the textural information, the larger the entropy value. The second moment was used to describe the distribution of image brightness. Correlation reflects the interrelation between the target pixel and its neighbors, and the correlation is strong if the interdependence is large.

3.2. Construction of Environmental Variables

The crop growth characteristics significantly vary under different environmental conditions. In the study area, the main environmental variables [29,30] affecting crop growth were as follows: (1) temperature, (2) precipitation, (3) slope, and (4) aspect (see

Table 3. Environmental indicators highly related to crop growth.

Index	Meaning
Slope	Degree of surface inclination
Aspect	Orientation of the topographic slope
$Y-Precipitation$	Mean precipitation during phenological period
$Y-Temperature$	Mean temperature during phenological period

). Temperature and precipitation were downloaded from WorldClim 2 [21]. Slope and aspect were calculated on the basis of the Shuttle Radar Topography Mission [31]. Details of data processing and calculation of environmental variables were obtained from previous studies [32,33].

3.3. ACCSH Method

This study proposed an ACCSH, which consists of three separate phases, as described below in Section 3.3.1, Section 3.3.2 and Section 3.3.3. Phase 1 aims to mine the spatial heterogeneity patterns of crop growth characteristics under the effects of multiple environmental variables. An adaptive spatial heterogeneity mining method (SHGD) presented in Section 3.3.1 was put forward on the basis of the geographical detector technique. In Phase 2, the classification indexes were optimized. As the crop growth characteristics differ in each spatial homogeneous zone obtained via SHGD, the key indexes for classification vary from zone to zone. Therefore, a separate index optimization in each spatial homogeneous zone must be performed to obtain effective classification indexes and reduce data redundancy. The optimization strategy is described in Section 3.3.2. In Phase 3, we constructed the classification model while considering the spatial heterogeneity, as explained in Section 3.3.3. The procedure is given in Figure 5.

3.3.1. Spatial Heterogeneity Patterns Mining

The geographical detector technique was proven effective for mining the spatial heterogeneity patterns under the effects of driving variables [10,11], but only discrete variables were suitable. Additionally, with numerous driving variables and categories, the spatial homogeneous zones may be too broken, which is inconvenient for subsequent classification operations. We addressed the abovementioned problems using an improved SHGD. The flow chart of SHGD is given in Figure 6, and its detailed procedure is as follows:

Step 1: Discretize continuous variables.

Environment variables driving the spatial heterogeneity of crops are mainly continuous measurements. The K-means method can segment continuous variables into several classes with similar intra-class distributions and different inter-class distributions. However, the number of classes in this method needed to be set. The NbClust package, which was developed to choose the optimal number of clusters for the dataset, was verified to be effective [34,35]. Hence, this study used the K-means function in the NbClust R package to obtain the optimal discrete categories.

Step 2: Determine key driving environmental variables.

The factor detector of the GD was adopted to mine the key driving environmental variables. However, highly correlated driving environmental variables need optimization. The detailed operations are as indicated below.

Step 2.1: In Equation (1), environmental variables are labeled as $X=\{x_1,x_2,\cdots ,x_m\}$ and the main classification index measuring crop growth characteristics is set as $Y$. Temporal NDVI can quantify crop growth characteristics to a large extent [36,37], and therefore, the mean NDVI of crop phenology is calculated as $Y$ in Equation (1). The environmental variables with a significant p-value less than 0.05 of q in Equation (1) are selected as key factors using the GD [38].

$$\left\{\begin{array}{c}q\left(x_i\right)=1-\frac{SSW}{SST}i=1,\cdots ,m\\ SSW=\sum _{h=1}^L{N}_h{\sigma }_h^2\\ SST=N{\sigma }^2\end{array}\right.,$$

(1)

where $m$ is the number of potential driving environmental variables $X$; $L$ is the class number of the environment variable $x_i$; $N_h$ is the number of samples in class $h$ of $x_i$; ${\sigma }_h^2$ are the variance of $Y$ in the area of class $h$ of $x_i$; ${\sigma }^2$ is the variance of $Y$ in the whole region; and $SSW$ and $SST$ are the within and the total sum of squares, respectively. The value range of $q$ is [0, 1], and a large $q$ indicates the strong ability to reveal the spatial heterogeneity of $Y$.

Step 2.2: Environmental variables are optimized. The Pearson’s correlation coefficients between key factors obtained in step 2.1 were calculated. Two variables with a Pearson’s correlation coefficient larger than 0.8 and a significance value smaller than 0.05 were regarded as highly correlated [39]. If two variables are highly correlated, only the variable with higher q-statistic significance in the GD result is retained. The reserved variables are labeled as $V$.

Step 3: Merge classes of environmental variables to eliminate invalid division.

If classes $A$ and $B$ simultaneously meet the following rule, they can be merged.

Rule: Classes $A$ and $B$ are from the same environmental variable. $Y$ in class $A$ is similar to that in class $B$, which can be measured using a t-test value of the risk detector of the GD result.

Step 4: Interaction of environmental variables. The detailed procedure is as follows.

Step 4.1: Classes in $V=\{v_1,v_2,\cdots v_k\}$ are labeled as $v_1=\left\{C_{v_1-1},C_{v_1-2},\cdots C_{v_1-n_1}\right\},v_2=\left\{C_{v_2-1},C_{v_2-2},\cdots C_{v_2-n_2}\right\},\cdots ,v_k=\{C_{v_k-1},C_{v_k-2},\cdots C_{v_k-n_k}\}$.

Step 4.2: Select two variables $v_i$ and $v_j$ to construct interactive classes. $v_i=\{C_{v_i-1},C_{v_i-2},\cdots C_{v_i-n_i}\}$ and $v_j=\{C_{v_j-1},C_{v_j-2},\cdots C_{v_j-n_j}\}$ can be interacted with as new classes $v_{ij}=\{C_{v_i-1\cap v_j-1},C_{v_i-2\cap v_j-1},\cdots C_{v_i-n_i\cap v_j-n_j}\}$. $v_{ij}$ is placed into $V$, and $v_i$ and $v_j$ are deleted from $V$.

Step 4.3: Go to step 4.2. The procedure is continued until only one element remains in $V$.

Step 5: Merge classes obtained by step 4. The difference of $Y$ between classes in $V$ is assessed using the t-test value of the risk detector of the GD result. Classes with no significant difference of $Y$ are merged.

Step 6: Calculate the spatial homogeneous zones. Locations of classes that remain after step 5 are the eventual spatial homogeneous zones $S=\{s_1,s_2,\cdots ,s_k\}$, where $k$ is their number. In these zones, $Y$ have remarkable differences driven by multiple environmental variables. That is, crop growth characteristics significantly differ between zones, and the samples were extracted for the following independent classification operation.

3.3.2. Index Optimization

Index optimization is generally carried out from two aspects. The first is to reduce the redundancy of indexes, commonly performed by assessing Pearson’s correlation coefficient. The second is to retain importance indexes to improve the efficiency of the classification model. RFAA is one of the typical methods for retaining importance indexes. RFAA can provide a feature ranking based on the contribution of each feature and selects optimal indexes by measuring the change in model performance when the order of each feature is disrupted. RFAA has been widely used and proven effective for index optimization [40,41].

In this study, correlations may exist among the spectral and textural indexes of the crops and indicate that indexes contain duplicate information or are prone to redundancy. In addition, the growth characteristics between a particular crop and other types are concentrated on specific spectral or textural features in certain months. Therefore, to further effectively reduce redundancy and improve accuracy, we optimized the classification indexes using not only Pearson’s correlation coefficient but also RFAA. Additionally, given the distinct variabilities in crop growth characteristics within each spatial homogeneous zone, index optimization was carried out individually. The detailed procedure of the index optimization strategy is explained below (Figure 7).

Step 1: Select one unlabeled spatial homogeneous zone for labeling.

Step 2: Perform correlation screening of indexes.

A Pearson’s correlation analysis was carried out for spectral and textural indexes. The correlation between indexes was considered high when the coefficient was larger than 0.8, and the significance value was smaller than 0.05 [39]. When indexes were highly correlated, only one index was retained.

Step 3: Perform optimal classification index mining.

RFAA is widely used in index screening, in which indexes are selected by measuring the variation of model performance when the order of each index is disrupted. Indexes were selected by evaluating their importance. The optimal combination of indexes can contribute the most to crop classification and was retained [42,43].

Step 4: Go to step 1 until classification indexes in all spatial homogeneous zones are optimized. Optimized indexes in each spatial homogeneous zone were saved and used for the following operations.

3.3.3. Classification Model Construction

The crop growth characteristics differ between spatial homogeneous zones, and the main indexes distinguishing between crop types vary in each spatial homogeneous zone. Hence, the crop classification must be separately carried out in each spatial homogeneous zone.

RF is a comprehensive learning algorithm based on decision trees, which has been widely used in various classification and regression fields [44,45]. The RF classification model was constructed using the training set in each spatial homogeneous zone based on the optimized indexes. The classification accuracy was evaluated using the validation set in each separate spatial homogeneous zone. The validation set construction strategy and classification evaluation indexes are given in Section 3.4.

3.4. Accuracy Evaluation

The Kennard–Stone (KS) algorithm can effectively select the training and validation sets with sufficient regional representation of classification indexes and was adopted in this study [46]. Among the sample crops, 70% were chosen as the training set, and the remaining 30% were utilized as the validation set [47]. The classification performance was eventually evaluated using the kappa coefficient, F1 score, and accuracy, which are commonly used and verified to be effective [48]. The F1 score is aimed at dichotomies and reflects the classification effect of target-type crops, which is the harmonic value of recall and precision. Additionally, accuracy signifies the percentage of the correctly predicted rate of the whole sample [49]. The kappa coefficient is a comprehensive index to evaluate the classification performance of the model as a whole.

4. Results

Wheat and corn are two crop types used in this study, and the classification results using ACCSH are given in Section 4.1 and Section 4.2, respectively.

4.1. Crop Classification of Wheat

The SHGD in ACCSH is used to mine the spatial heterogeneity patterns of wheat. The mean NDVI is used as Y in Equation (1). Factor detection using steps 1 and 2 of SHGD indicates that $Y-Precipitation$, which is divided into three classes, has a significant effect on the spatial heterogeneity characteristics of wheat growth (

Table 4. Factor detection result using the GD model for wheat classification.

Statistical Variable	$\mathit{Y}-\mathit{P}\mathit{r}\mathit{e}\mathit{c}\mathit{i}\mathit{p}\mathit{i}\mathit{t}\mathit{a}\mathit{t}\mathit{i}\mathit{o}\mathit{n}$	$\mathit{Y}-\mathit{T}\mathit{e}\mathit{m}\mathit{p}\mathit{e}\mathit{r}\mathit{a}\mathit{t}\mathit{u}\mathit{r}\mathit{e}$	Aspect	Slope
q statistic	0.10	0.01	0.01	0.00
p-value	0.00	0.14	0.43	0.94

). Risk detection using step 3 of SHGD in

Table 5. Risk detection result of classes of $Y-Precipitation$ using the GD model for wheat classification.

Classes	Class 1	Class 2	Class 3
class 1
class 2	Y
class 3	Y	N

Y represents the remarkably different crop growth characteristics between classes and vice versa anti.

indicates that the crop growth characteristics of wheat in classes 2 and 3 of $Y-Precipitation$ have no significant difference. As such, classes 2 and 3 of $Y-Precipitation$ can be combined as shown in

Table 6. Optimized classes of $Y-Temperature$ for wheat classification.

Classes	Optimized Classes	Values	Spatial Homogeneous Zones
class 1	A	$Y-Precipitation$ between 43.0 and 57.0	I
class 2	B	$Y-Precipitation$ between 57.0 and 68.8	II
class 3	B	$Y-Precipitation$ between 68.8 and 109.6	II

, and the whole area is divided into spatial homogeneous zones I and II, respectively (Table 6 and Figure 8).

According to Section 3.1, spectral and textural indexes during the phenological period (Figure 3) are calculated for wheat classification. The index optimization strategy described in Section 3.3.2 is carried out on the spectral and textural indexes. The results in

Table 7. Optimized classification indexes of wheat.

Region	Index Type	Index Name
Global region	Spectral indexes	1-NDVI, 5-NDVI, 7-NDVI, 8-NDVI
		5-BLUE, 8-BLUE
		7-NIR,12-NIR
		8-MIR
		$T_i$, $D_i$
	Textural indexes	6-B4, 7-B4, 12-B4
		12-B5
		7-B6
		6-B7
Zone I	Spectral indexes	3-NDVI, 4-NDVI, 5-NDVI, 7-NDVI, 8-NDVI
		10-RED
		4-NIR
		5-MIR
		$T_i$, $D_i$
	Textural indexes	8-B1
		9-B2
		10-B6
		9-B8, 12-B8
Zone II	Spectral indexes	5-NDVI, 7-NDVI, 8-NDVI
		7-RED, 8-RED, 9-RED
		8-BLUE, 9-BLUE
		8-MIR, 9-MIR
		$T_i$, $D_i$
	Textural indexes	7-B1, 8-B1
		7-B4, 8-B4
		7-B6, 8-B6
		8-B7, 9-B7

show that the optimal index sets of wheat vary within each spatial homogeneous zone, further verifying the necessity of index optimization.

The ACCSH based on the optimized classification indexes is carried out, and the corresponding classification evaluation indicators are shown in Figure 9. The ACCSH in zones I and II show F1 scores of 0.94 and 0.95, respectively; the values of classification accuracy are 0.93 and 0.93, respectively, and the kappa coefficients are 0.86 and 0.84, respectively. The evaluation indicators indicate that ACCSH exhibits high stability and classification ability and can be utilized to classify wheat in a large-scale region with significant spatial heterogeneity.

4.2. Crop Classification of Corn

SHGD is used to mine the spatial heterogeneity patterns, and the mean NDVI of corn during the phenological period is used as $Y$ in Equation (1). The results indicate that $Y-Precipitation$ and $Y-Temperature$ can be divided into three classes using K-means. Then, key driving environmental variables are obtained, and both $Y-Precipitation$ and $Y-Temperature$ show significant effects (p < 0.05) on the spatial heterogeneity of corn growth characteristics according to the factor detection result in

Table 8. Factor detection result using the GD model for corn classification.

Statistical Variable	$\mathit{Y}-\mathit{P}\mathit{r}\mathit{e}\mathit{c}\mathit{i}\mathit{p}\mathit{i}\mathit{t}\mathit{a}\mathit{t}\mathit{i}\mathit{o}\mathit{n}$	$\mathit{Y}-\mathit{T}\mathit{e}\mathit{m}\mathit{p}\mathit{e}\mathit{r}\mathit{a}\mathit{t}\mathit{u}\mathit{r}\mathit{e}$	Aspect	Slope
q statistic	0.01	0.04	0.02	0.00
p-value	0.02	0.02	1.00	0.83

. According to step 3 of SHGD, the risk detection results are given in

Table 9. Risk detection result of classes of $Y-Precipitation$ using the GD model for corn classification.

Classes	Class 1	Class 2	Class 3
class 1
class 2	Y
class 3	Y	N

and

Table 10. Risk detection result of classes of $Y-Temperature$ using the GD model for corn classification.

Classes	Class 1	Class 2	Class 3
class 1
class 2	N
class 3	N	Y

, respectively, and classes of $Y-Precipitation$ and $Y-Temperature$ are optimized as classes in

Table 11. Optimized classes of $Y-Precipitation$ for corn classification.

Classes	Optimized Classes	Values
1	C	$Y-Precipitation$ between 48.5 and 60
2	D	$Y-Precipitation$ between 60 and 71.8
3	D	$Y-Precipitation$ between 71.8 and 95.5

and

Table 12. Optimized classes of $Y-Temperature$ for corn classification.

Classes	Optimized Classes	Values
1	E	$Y-Temperature$between 10.0 and 13.3
2	E	$Y-Temperature$ between 13.3 and 14.2
3	F	$Y-Temperature$ between 14.2 and 15.5

, respectively, and are further interacted according to step 4. The risk detection of interaction classes is subsequently executed (

Table 13. Risk detection result of interaction classes using the GD model for corn classification.

Classes	Class 1	Class 2	Class 3	Class 4
class 1
class 2	Y
class 3	Y	Y
class 4	Y	Y	Y

), which manifests their remarkable variation in crop growth characteristics. Therefore, the locations of interaction classes are considered as the eventual spatial homogeneous zones (

Table 14. Interaction classes of $Y-Precipitation$ and $Y-Temperature$ for corn classification.

Classes	Interaction Information	Values	Spatial Homogeneous Zones
1	$\mathrm{C}\cap \mathrm{E}$	$Y-Precipitation$ between 48.5 to 60 and $Y-Temperature$between 10.0 and 13.3	Zone I
2	$\mathrm{C}\cap \mathrm{F}$	$Y-Precipitation$ between 48.5 to 60 and $Y-Temperature$ between 14.2 and 15.5	Zone II
3	$\mathrm{D}\cap \mathrm{E}$	$Y-Precipitation$ between 60 to 95.5 and $Y-Temperature$between 10.0 and 14.2	Zone III
4	$\mathrm{D}\cap \mathrm{F}$	$Y-Precipitation$ between 60 to 95.5 and $Y-Temperature$between 14.2 and 15.5	Zone IV

and Figure 10).

Spectral and textural indexes of corn during its phenological period are calculated according to methods in Section 3.1. Similar to the results for wheat, the optimized classification indexes in spatial homogeneous zones and global areas vary (

Table 15. Optimized classification indexes of corn.

Region	Index Type	Index Name
Global region	Spectral indexes	5-NDVI, 7-NDVI
		7-RED, 8-RED
		5-BLUE
		9-NIR
		5-MIR
		$T_i$, $D_i$
	Textural indexes	5-B1, 6-B1, 12-B1
		7-B3
		7-B7
		6-B8
Zone I	Spectral indexes	5-NDVI, 7-NDVI, 8-NDVI
		5-RED, 8-RED
		5-MIR, 6-MIR
		$D_i$
	Textural indexes	1-NDVI, 2-NDVI, 3-NDVI, 7-NDVI, 11-NDVI, 12-NDVI
		2-B6
		7-B7
		4-B8, 6-B8
Zone II	Spectral indexes	5-NDVI, 6-NDVI, 7-NDVI, 8-NDVI
		7-RED
		7-BLUE, 8-BLUE
		7-NIR
		7-MIR
	Textural indexes	2-B1, 5-B1, 6-B1, 12-B1
		1-B2
		2-B4
		1-B6
		1-B7
Zone III	Spectral indexes	11-NDVI
		5-RED, 7-RED, 8-RED
		4-BLUE, 5-BLUE, 8-BLUE
		7-NIR
		4-MIR, 5-MIR, 8-,MIR
		$D_i$
	Textural indexes	5-B1
		2-B2
Zone IV	Spectral indexes	6-NDVI, 8-NDVI
		4-RED, 5-RED, 6-RED, 7-RED, 8-RED, 11-RED
		4-BLUE, 5-BLUE
		7-NIR
		5-MIR
		$D_i$
	Textural indexes	3-B1, 7-B1
		7-B3

), setting barriers to the high accuracy of crop classification at the global scale.

Figure 11 presents the classification evaluation indicators of corn. Results show that the F1 scores and accuracy indicators of ACCSH in zones I–IV are all higher than 0.9, verifying its effective crop classification ability with high accuracy in a large-scale region with significant spatial heterogeneity.

5. Discussion

5.1. Comparison Results of ACCSH with Classical Methods for Crop Classification

The classification effectiveness of ACCSH’s zoning strategy is further verified by comparison at the global scale; that is, the classification is constructed using RF combined with optimized classification indexes at the global scale in the study area. The validity of the index optimization operation can be verified using ACCSH and ACCSH with non-optimized indexes. Additionally, one of the widely used deep learning methods, the 1D Convolutional Neural Network (1D-CNN), is used for comparison. The comparative results are given and analyzed below.

Figure 12a shows the evaluation indicators of wheat. The results show that the kappa coefficients of ACCSH in zones I and II are 16% and 14% higher than that of the classification at the global scale, respectively. Figure 12b presents the classification evaluation indicators of corn. The results show that the kappa coefficients of ACCSH in zones I–IV are all higher than those at the global scale. Thus, Figure 12 indicates that the index optimization strategy combining correlation screening and RFAA in Section 3.3.2 can effectively improve crop classification accuracy. Classification based on an index set optimized using correlation screening can improve the accuracy compared with that based on a non-optimized index set. Accuracy can be further improved based on the index set optimized using correlation screening followed by RFAA, and the index optimization strategy is adopted in ACCSH as described in Section 3.3.2. The index optimization strategy adopted in ACCSH can improve the kappa coefficient of classification effect by 5–18% for wheat and 9–42% for corn. The combination of correlation screening and the RFAA index optimization strategy can effectively screen the key classification indexes. In summary, the index optimization and spatial zoning strategies of ACCSH can improve crop classification. The kappa coefficient of ACCSH is improved by 20–67%, and the F1 score of ACCSH is improved by 6% compared with that of the global classification and classification based on non-optimized index sets, respectively. The ACCSH can perform crop classification with relatively high ability and can meet actual demand in large areas with significant spatial heterogeneity.

The kappa coefficients of wheat and corn using 1D-CNN classification are 0.72 and 0.58, respectively, which are significantly lower than those of ACCSH. The possible reason is that the training of deep learning models requires a large number of datasets [50,51]. When the sample size is limited, the machine learning-based ACCSH can achieve a relatively good classification accuracy.

5.2. Analysis of the Spatial Heterogeneity Patterns of Crop Growth

According to Section 4.1, precipitation has a significant effect on wheat growth. The mean annual precipitation in zone II is higher than that in zone I (Table 4), and the mean NDVI value of wheat in zone II is higher than that in zone I during the planting and growth period of wheat (Figure 13a). That is, in the study area, precipitation has a positive driving effect on wheat growth, which may be the reason why a certain level of high precipitation results in better soil moisture reserves, contributing to better conditions. According to Section 4.2, the growth characteristics of corn significantly vary within spatial homogeneous zones driven by precipitation and temperature (Figure 13b). The higher precipitation during the growth period of corn in zones III and IV causes better growth than in zones I and II, and this finding is highly consistent with the relationship between precipitation and wheat growth. Additionally, for corn, the higher temperature in zones II and IV reduces the NDVI values compared with zones I and III during the planting period. High temperature increases the evaporation of water and reduces humidity, resulting in difficulties for emerging seedlings and the water loss of pollen. The interaction of these two factors shows that the combination of high precipitation and low temperature in zone III contributes to the best corn growth, whereas the combination of low precipitation and high temperature are the worst conditions. In summary, SHGD in ACCSH can mine the spatial heterogeneous characteristics of crop growth. To further demonstrate the efficacy of SHGD, we use the typical layered approaches, such as that based on management zones [7,8] and on segmented single climatic variables (i.e., temperature or precipitation) in place of SHGD in the ACCSH. The results in Figure 14 indicate that ACCSH based on SHGD has the highest classification accuracy, demonstrating its power to mine spatial heterogeneous patterns for crop classification. Specifically, the best wheat zoning outcome depends on the precipitation segmentation that SHGD adaptively mines.

In summary, ACCSH can be effectively applied to crop classification in large areas, and the proposed SHGD can mine the driving effects of multiple environmental variables on the spatial heterogeneity of crop growth characteristics.

6. Conclusions

This study proposes ACCSH to automatically classify crops considering the spatial heterogeneity of crop growth characteristics under the influence of multiple environmental factors. Experiments on two real applications are carried out. Comprehensive comparisons are designed, including classifications at the global scale and based on non-optimized classification indexes. The proposed ACCSH has the following advantages: (1) The results of SHGD in ACCSH can reveal the spatial homogeneous zones and mechanism of crop growth characteristics under the influence of multiple environmental variables. (2) The index optimization strategy in ACCSH can fully consider the differences between spatial homogeneous zones by carrying out a separate optimization in each spatial homogeneous zone and obtaining the key classification indexes to improve classification accuracy. (3) ACCSH can automatically perform crop classification without sufficient prior knowledge of the study area.

ACCSH is used for wheat and corn classification in northern France. The following new findings are obtained from the study: (1) The influence of the mean annual precipitation on wheat is significantly positive. The better condition for wheat growth is due to relatively high precipitation; (2) Corn growth is sensitive to both precipitation and temperature, where the former exhibits an overall positive driving effect, whereas the latter tends to have a negative effect on crop growth; (3) The accuracies of ACCSH evaluated using the kappa coefficient are 15% for wheat and 26% for corn, higher than those of classification at the global scale; (4) The accuracies of ACCSH are also much higher than those of classifications based on non-optimized classification indexes.

Overall, ACCSH can automatically carry out crop classification under the influence of multiple environmental factors on crop growth, indicating its potential in crop classification with relatively high precision. In addition, ACCSH is adaptive and implemented using R 4.3.0 software and Matlab 2014a software, which makes the operation convenient. Future research can focus on the application and improvement of the ACCSH. The proposed method is suitable for single crop classification, and future studies can improve ACCSH to adapt to multiple and rotation crop classifications. Additionally, if sufficient remote sensing data with high spatial resolution are available, ACCSH can be carried out at a small scale.