Crop classification based on the spectrotemporal signature derived from vegetation indices and accumulated temperature

ABSTRACT Due to differences in environmental factors, the phenology of the same crop is different every year, causing divergent performances of the classifier built by spectral or time-series features Here, we proposed a random forest classifier (RFC) based on an asymmetric double S curve model fitted by accumulated temperature (AT) and Vegetation Index (VI), which can be applied in different years without ground samples. We built AT and VI time series from Moderate Resolution Imaging Spectroradiometer 8-day composites of land surface temperatures and Sentinel-2 and Landsat-8, respectively. The RFC was trained by characteristics from the asymmetric double S curve. We prepared RFC by ground samples of 2018 and 2019 and then mapped crops of the same region in 2017. Results indicated that, compared with diverse VI-AT series, the overall accuracy based on universal normalized vegetation index (UNVI) was the best of all (2017: F1 = 0.91, 2018: F1 = 0.92, 2019: F1 = 0.91) and better than that based on the UNVI-TIME series (2017: F1 = 0.84, 2018: F1 = 0.81, 2019: F1 = 0.88). It proved that the classification features from the VI-AT series have smaller intra-class differences in 2017, 2018, and 2019.


Introduction
As one of the most important corn and soybean producers, the USA produces one-third of the world's corn production. Most of the corn is planted in a belt region, which extends across 12 Midwest states and is called Corn Belt (Panagopoulos et al. 2015). Due to the influences of climate, cost of inputs, cropping patterns, and other factors, the cropping area changes every year, alternating planting corn and soybean. The corn cropping area could extend to the northwest of the USA (Green et al. 2018). This could lead to historic carbon losses, and no-till practices could help restore the lost soil organic carbon (Mishra, Ussiri, and Lal 2010;d'Andrimont et al. 2020). Insecticide applications may cause environmental pollution (Hladik, Kolpin, and Kuivila 2014). Therefore, it is important to develop a method to detect corn area change for agricultural management and environmental protection.
Over the past 20 years, significant progress has been made in the application of remote sensing, and this technique has been widely used in crop classification, cropping area monitoring, and crop yield prediction (Zhao 2014;Wójtowicz, Wójtowicz, and Piekarczyk 2016). Crop classification is the base of most remote-sensing applications in agriculture.
The classifier based on a single-date image was first used for land cover classification. But the different landcover with the same spectrum at a given period limits the classification accuracy. To overcome this problem, classification based on time-series remote sensing has been proposed. Compared with single-date remote sensing images, time-series remote sensing contains more information about a crop's biophysical features derived from multi-temporal images. In particular, time series can reveal spectral changes reflecting changes in phenology (Bai et al. 2020). Time series have been widely used in crop classification, and the accuracy of the classifications based on time series is better than that of classifications based on mono-temporal remote sensing (Turker and Arikan 2005;Ozdogan and Gutman 2008;Hao, Niu, and Wang 2014;Liu et al. 2020b;). In the following research, time-series data played an important role in agriculture and ecology. Krehbiel, Jackson, and Henebry (2016) monitored the influence of the urban heat island effect on land surface phenology. Wang et al. (2020b) mapped the area of sugarcane using time series data derived from Sentinel-1, Sentinel-2 and Landsat. Shen et al. (2019) used a 10-year time series to identify the blooming of various seaweeds. Chen et al. (2020a) monitored changes in the winter wheat planting area via time-series remote sensing. Early crop monitoring based on time series also has been intensely investigated (Johnson and Mueller 2021;Hao, Di, and Guo 2022).
Improvements in the classification method, such as the development of machine learning (ML) (Wang et al. 2020a;Yu et al. 2020) and deep learning (DL) (Li et al. 2017Liu et al. 2020a;Liu, Yang, and Lunga 2021;Xu et al. 2021), have enabled to capture sufficient characteristics from timeseries images. Bofana et al. (2020) compared crop classifications using the four main ML methods based on time series and analyzed the benefits of different methods. The classification based on ML and DL has been used in agricultural mapping and monitoring (Kussul et al. 2015;Xu et al. 2020) and land monitoring (Roberto et al. 2019). Although classification based on ML and DL has higher classification accuracy than traditional classification, there are many challenges while applying these methods; the transformation of classifiers limited the applicability of ML and DL. In most cases, there are insufficient valid training samples for a given time, and the performance of classifiers degenerates when it is used for a new time period.
For solving the problem, we used the relation between Land Surface Temperature (LST) and Vegetation Index (VI) as the classification feature to build a robust classifier in different years. Some studies have attempted to classify crops by combining environmental factors with spectral or time-series characteristics to enhance the classifier performance (Foerster et al. 2012). de Beurs and Henebry (2004) used the accumulated growing degree-day (AGDD) to associate temperature with the normalized difference vegetable index (NDVI). Nguyen et al. (2020) monitored the transformation between grass and farmland in the Corn Belt by incorporating the AGDD and NDVI. Hou et al. (2014) used the growing degree day (GDD) to analyze the influence of temperature on crop growth. Qian et al. (2019) studied the relationship between the growing degree unit (GDU) and crop growth. Compared with AGDD, GDD, and GDU have different conditions, limiting the temperature interval to 0-30°C during the calculations of AT. Eisavi et al. (2015) found that the combination of temperature and time-series data can improve classification accuracy. Chen et al. (2020b) found that combining modified LST and NDVI time series enhanced the classifier's performance. However, it has been effective in improving classification accuracy by combining the environment factor and spectral feature. On the one hand, the environmental factor was often used as an assistant in part of studies; on the other hand, the relation between the LST and VI was shown by a quadratic function, and it is hard to reflect the characteristic of different crops.
In this study, we used the crop area in the USA Corn Belt as the study area and built the asymmetric double-sigmoid function through VI time series and LST time series. In this way, we extracted the stable crop characteristics in different years and built a crop classifier that can keep classification accuracy in different years without transforming training. Our work provides a new method to improve classification applicability in different years and quantify the response of vegetation to the change of environment. Meanwhile, our study offered basic research for early crop detection and crop health analysis based on environmental factors in the future.

Study area
The areas of Iowa and Nebraska mixed with corn and soybean in the USA are chosen as the study areas for this research (Figure 1), because they are traditional agricultural states, and corn and soybean are the main agricultural outputs. They have a continental climate with considerable seasonal variation, which was likely to enhance the separability of different crops based on the use of the temperature field. In addition to farmland, there is also non-agricultural land cover, such as buildings, grass, and forest in the study area.
The area in Iowa was approximately 21 km wide and 38 km long, and the area in Nebraska was approximately 28 km wide and 39 km long. According to the United States Department of Agriculture (USDA) Cropland Data Layer (CDL) (Boryan et al. 2011), approximately 80% of the farmland was used for corn and soybean in 2016 in both study areas, and the remainder comprising small areas of alfalfa, wheat, canola and wetland, and wetland, which were located along the river. According to the main crop calendar from USDA (Figure 2), there are conspicuous differences between the main crop phenology. Significantly, the corn planting period is longer than the soybean planting period. Meantime, the crop calendars in 2019 are lower than those in 2018. We will research the classification based on the data of the Iowa study area in 2017, 2018, and 2019 and verify the classification method based on the data of the Nebraska study area in 2017, 2018, and 2019.

Data and pre-processing
We built NDVI, UNVI, and NDWI time series derived from Sentinel-2 and Landsat-8 data, and the accumulated temperature (AT) time series based on Moderate Resolution Imaging Spectroradiometer (MODIS) data.

Vegetation Indices (VI)
Various VIs can reflect the change of vegetation phenology by the relation between reflectivity and physical characteristics of vegetation (Jiang et al. 2008;Brown et al. 2013;Gong et al. 2013;Testa et al. 2018;Seo et al. 2019).
As the representation of all VI, NDVI has been widely used after its establishment in 1973 (Rouse 1973). The NDVI was calculated using Landsat-8 and Sentinel-2 data as shown below: where r NIR is the reflectance in the NIR band of Landsat-8 and Sentinel-2 and r RED is the reflectance in the red band of Landsat-8 and Sentinel-2. Compared with VIs calculated directly by several separate bands, UNVI was calculated by the parameters extracted by the whole bands to make the best use of the spectral information (Zhang et al. 2007). For the UNVI, the surface reflection of a pixel was referred to as a linear combination of water, vegetation, soil, and yellow-leaf parameters, as follows: where i represents the band number of a given sensor, R i represents the pixel reflection in the band i, P w(i) , P v(i) , P s(i) , P 4(i) are the standard spectral patterns of water, vegetation, and soil and the supplementary yellow-leaf pattern in the band i, respectively; and C w , C v , C s and C 4 are the decomposition coefficients in the band i. The last function can be simplified as shown below: where C = [C w, C v , C s , C 4 ], M is the constant matrix determined by the sensor type, and R is the spectral data of every band in the pixel. For Multi-Spectral Imager (MSI) on Sentinel-2 and Operational Land Imager (OLI) on Landsat-8, M is determined as follows: The UNVI can then be calculated using the following formula with parameters C w , C v , C s , and C 4 : In most cases, the constant term a is 0.1 . Apart from using UNVI and NDVI to show the vegetation growth, other VIs can reflect environmental condition of vegetation growth. The NDWI was proposed in 1996 to detect open-water features and study water quality (McFEETERS 1996). The NDWI is calculated as follows: where r GREEN is the reflectance in the green band and r NIR is the reflectance in the NIR band.

Pre-processing of the VI time series
To reduce the influence of the environment and ensure that the VI data corresponded to the LST data, we rebuilt the VI time series. In our study, we resized the time gap to 8-day, the same as the time resolution of LST data, by the following method.
where t is the time point of interest, V t is the value at time t, t b is the preceding time point, t a is the next time point after t, V t b is the value at time t b , and V t a is the value at time t a . In this way, the abnormal lower value in time series, led by cloud or shade, can be rebuilt by the value of adjacent time.
We obtained the final VIs time series by a Savitzky-Golay (SG) filter, which was widely used in previous studies and performs well (Chen et al. 2004;Tang et al. 2019). To maintain the change characteristics of the VI time series, we set the window of the SG filter to 5 and the polynomial order of the SG filter to 3 after comparing the version parameter setting.

Accumulated temperature time series
AT plays an important role in agricultural research and impacts vegetation growth. Due to the limitation of meteorological station distribution, most researchers used LST derived from remote sensing instead of the air temperature to calculate AT. The relationship between AT and the growth stage of vegetation is not linear. All vegetation has a base temperature span in which the vegetation can thrive in. The vegetation will stop growing if the temperature is above the maximum base temperature or below the minimum base temperature (Sacks and Kucharik 2011). Hence, we selected the Figure 3. Comparison between the NDVI time-series based on only Sentinel-2 data and that based on Sentinel-2 and Landsat-8 data for corn and soybean, respectively. The broken curves are NDVI time series based on Sentinel-2 and Landsat-8 data and the full curves are NDVI time series based on only Sentinel-2 data. The separability between soybean and corn NDVI time series based on Sentinel-2 and Landsat-8 data is better than that based on only Sentinel-2 data.
following method for calculating AT: where AT t is the AT at time t, AT t−1 is the AT at time t-1, T t is the temperature at time t, T t−1 is the temperature at time t-1, T min is the minimum base temperature, T max is the maximum base temperature. Based on the previous research (Hou et al. 2014;Qian et al. 2019), the T min was set as 0 and T max was set as 30 in our study. The Dt is the time gap in the time series. We used the Land Surface Temperature and Emissivity 8-Day Global 1-km LST product from MODIS to build the AT time series. The LST data were gained from Google Earth Engine (GEE), and the data were resampled to 10 m by the GEE resample method before downloading.
The LST can be calculated using the brightness temperature band of Landsat-8 or Landsat-7, which has a higher spatial resolution than MODIS. But limited by the number of the available images and the difference between the surface temperature derived from MODIS and Landsat (Figure 4), it is hard to calculate a fine AT time series based on the Landsat data.
Although the tendencies of the two surface temperature time series are similar, the land surface temperature derived from Landsat is greatly different from that derived from MODIS from day 200 to day 250. And it is tough to estimate truth temperature in this period based on Landsat (Figure 4).
Impacted by the rough space resolution of AT time series, the mixed pixels may influence the separability of different crops in the same year. But in our study, the separability of different crops is determined by the VI time series, derived from Sentinel-2 and Landsat-8, and their resolutions are 10 and 30 m, respectively. The difference of the same crop in different years is reduced using LST, and it can't be affected by the mixed pixels.

Ground reference
The CDL was produced by the USDA based on moderate resolution satellite imagery and extensive agricultural ground observations (Howard, Wylie, and Tieszen 2012;Hao et al. 2015;Li et al. 2015;Chen et al. 2020b;Momm, ElKadiri, and Porter 2020). It has been widely used in remote sensing as the truth ground data, and its accuracy of landcover can keep 80% in most areas. In our study, we randomly selected ground samples through study area CDL from 2017 to 2019.
In the study, 70% of each type of ground sample were selected randomly as a training sample, and the remaining ground samples were used as validation data for 2018 and 2019.
The term 'non-crop' was used to refer to vegetation cover except soybean and corn, mainly wetland, grassland, and forest. And other crops, such as alfalfa and oats, are too small to gain enough ground samples, so we also referred to them as non-crop. The statistical details of the ground samples in three years are shown in Table 3.

Relationship between phenology and temperature
In most studies, time is used as an indicator of the growth stage of vegetation (Sakamoto, Gitelson, and Arkebauer 2013;Vrieling et al. 2018;Konduri et al. 2020;Peng et al. 2021). However, the growth stage of vegetation is more strongly controlled by environmental factors than time. The relationship between vegetation and time is also controlled by the energy that vegetation absorbs from the environment over time. If environmental factors, such as surface temperature at the same time-point change, the growth stage will also change accordingly. Some studies have compared the changes in AT and changes in crop NDVI time-series (Qian et al. 2019) and shown that changes in surface temperature could impact plant phenology (Krehbiel, Zhang, and Henebry 2017). As shown in Figure 5, based on the comparison of the corn NDVI and AT time series between 2018 and 2019, it is evident that the faster the AT increases, the earlier the corn growth started. The surface temperature differences were consistent with differences in the NDVI time series.
To prove the efficiency of AT better, the NDVI time series of corn samples were divided into 28 sections from April 16 to November 18, while the AT time series was divided into 28 sections from 200°C to 5600°C in our study. The differences in the corn NDVI time series between 2018 and 2019 were then assessed by correlation coefficient. Figure 6 showed that there was little difference in the NDVI-AT series between 2018 and 2019. It is possible to improve the consistency of the NDVI time series in different years using AT as an indicator of vegetation growth.

The asymmetric double-sigmoid function
Many researchers have studied the relationship between surface temperature and vegetation growth and built the model to quantify the change of vegetation growth led by the surface temperature change (Dong et al. 2009;Baumann et al. 2017;Li et al. 2017;Hu et al. 2019;Bagherzadeh, Hoseini, and Totmaj 2020;Ji et al. 2021). In contrast, methods such as the double logistic model (Meroni et al. 2014) and the straightforward definition of the linear harmonic model (Beck et al. 2006) have been applied in previous studies. We used the asymmetric double-sigmoid function proposed by Zhong et al. (2011), because compared with the double logistic model or the straightforward definition of the linear harmonic model, the asymmetric double-sigmoid function would reveal spectral characteristics of vegetation during the nutrition-accumulation stage and provide a good balance between model complexity and the ability to reveal crop spectral characteristics. The  asymmetric double-sigmoid function is written as where N(t) is the VI value with the AT at time t, N b is the VI background value in the leafless season, N a is the amplitude of VI change in the growing cycle, D i is the AT at the maximum rate of increase in the VI, D d is the AT at the minimum rate of increase in the VI, and p and q are the rates of change of the increasing and decreasing slopes, respectively. Based on the VI time series and AT time series of every pixel, we used the non-linear leastsquares method to fit the asymmetric double-sigmoid function. Upon completion of the iteration, if the appropriate parameters for the pixel still could not be fitted, we set all parameters to 0.
To distinguish the curve based on fitted by different data efficiently, we called the fitted curve based on VI time series and AT time series as VI-AT series, and the fitted curve based on VI time series and TIME series as VI-TIME series.

Random forest classifier
The random forest classifier has been widely used in remote sensing classification (Zhang and Yang 2020;Dixon et al. 2021), and good performance has been reported compared with other ML methods, such as the support vector machine method . The random forest method has a higher classification accuracy for non-linear characteristics (Jhonnerie et al. 2015). In our study, the random forest classifier was used as a supervised classifier. The random forest classifier was trained based on fitted parameters (N b , N a , p, q and D d − D i ) from training points in 2018 and 2019.
The accuracy of the classification was assessed by the overall accuracy (OA) (Story and Congalton 1986;Congalton 1991;Naesset 1996;Stehman 1996;Rwanga and Ndambuki 2017) and F1-score (Fawcett 2006;Sasaki 2007). Figure 7 shows the flow chart used in the study. Firstly, VI time series were built using Sentinel-2 and Landsat-8 data, and AT time series were built using the Land Surface Temperature and Emissivity 8-Day Global 1-km product from MODIS. Next, we fitted the asymmetric double-sigmoid function through VI and AT time series. Finally, we built a classifier using parameters of the asymmetric double-sigmoid function and assessed the classifier performance by test points. To verify the classifier's performance in different years without training, we applied the classifier to map the crop field in the study area in 2017.

Mapping of crops
We also mapped the crop in the same way based on the VI-TIME series to study the influence of AT time series on the classification accuracy.

Classification results based on the VI-AT series for 2017, 2018, and 2019 in Iowa
We trained the random forest classifier by the sample points and VI-AT series of 2018 and 2019. We mapped the crop in 2017, 2018, and 2019. Table 4 contains the OA and Kappa coefficient of classification based on the NDVI-AT, UNVI-AT, and NDWI-AT series for 2017, 2018, and 2019.
Of the NDVI-AT, UNVI-AT, and NDWI-AT series, the classification based on the UNVI-AT series was most accurate in all periods. For the NDVI-AT and UNVI-AT series, there was little change in the classification accuracy between 2017, 2018, and 2019. It means there was a stable and similar crop characteristic based on the VIs-AT series in different years, and the characteristic based on the UNVI-AT series had better separability. The great change of classification accuracy based on the NDWI-AT series may be led by the error while fitting the Asymmetric Double-sigmoid Function.
As shown in Figure 8, the classification based on the NDVI-AT series exhibited the worst performance in classifying non-crop and city. In 2018, the city area was misclassified as non-crop cover, and the area of soybean fields was less than that in the CDL. In addition, some fields were mixed with soybean and corn in classification based on the NDVI-AT series in 2017; it is not the same as the ground reality. In 2019, parts of the non-crop area around watercourses were incorrectly classified as city areas. The classification based on the NDWI-AT series was the least accurate due to the bad fitting result in 2019. Some non-crop areas were incorrectly classified as city areas in 2018 and 2019, and the city area was largely misclassified.
Incorporating the UNVI time series and AT time series realized the advantages of a VI-AT series better than the NDVI and NDWI time series. In future studies, the UNVI-AT series will be used as the representative VI-AT series.

Classification results based on the UNVI time-series for 2017, 2018, and 2019 in Iowa
To highlight the advantages of the VI-AT series, the study also mapped the crop based on the UNVI-TIME series. The method used to train the classifier based on the UNVI-TIME series was the same as that of the UNVI-AT series. The accuracy of the classification is shown in Table 5.   The OA of classification based on the UNVI-AT series was better than that of the UNVI-TIME series for 2017, 2018, and 2019. There was little difference in the classification accuracy between the UNVI-AT series between 2017, 2018, and 2019. In contrast, there were relatively large differences in the classification accuracy based on the UNVI-TIME series between 2017, 2018, and 2019. Therefore, the classifier characteristics for the same kinds of land cover based on the UNVI-AT series were more similar than when the UNVI-TIME series was used. Figure 9 proved the classifier based on the UNVI-TIME series has worse performance in classifying city and non-corp. The area of the city in the classification based on the UNVI-TIME series was larger than the truth. The parts of the soybean area were classified as corn and city in 2018. And in 2017, the lower accuracy of classification was due to parts of the soybean field being misclassified as corn. In other words, the vegetation spectral characteristics were impacted by the change of environment.

The classification applied in Nebraska in 2017, 2018, and 2019
In our study, the classification method was also applied in the study area of Nebraska. The classification ( Figure 10) proved that the method could well classify most of the cornfields and soybean fields in 2017, 2018, and 2019. Parts of soybean or corn field were misclassified as alfalfa, and the accuracy of grassland and forest was not good.
The overall accuracy and Kappa of classification are shown in Table 6. According to the comparison between classification output and CDL, the accuracy change may be led by misclassification in terms of forest, grassland, and alfalfa.

Differences in classification accuracy based on the VI-AT series
In our study, we compared the classification based on various VI-AT series to study the characteristic change in different years. As shown in Table 7, the classification accuracy based on the UNVI-AT series was better than the other classification accuracy, not only for the overall accuracy but also for the single landcover.
According to the crop map based on the NDWI-AT series and the classification accuracy, there are several reasons for the poor performance of the NDWI-AT series. Because we set the initial parameters and number of iterations, and the pixel parameter was 0 if the good fitting parameters were not obtained after the iterations. From the fitting result based on the NDWI-AT series, there were many pixels in corn and soybean fields whose fitted parameters were 0 in 2019. This led many cornfields to be misclassified as soybean fields. The variance of the parameters fitted by the NDWI-AT series was larger than those fitted by the UNVI-AT and NDVI-AT series. It was, therefore, not possible to obtain the right parameters for most of the NDWI-AT series pixels with only 3000 iterations. The NDWI is typically used to show the water condition in a canopy rather than the vegetation growth state. The NDWI change may be related to the crop growth stage and the environmental factor, such as irrigation. The result from curve fitting proved that the relationship between the NDWI and AT was not the same as that between the UNVI and AT or the NDVI and AT.   As for the classification accuracy difference in the UNVI-AT series and NDVI-AT series, the advantage of the classifier using the UNVI-AT series is the method used to calculate the UNVI. The UNVI is calculated using all bandsfrom visible to short wave radiationand includes information about reflectance from water, soil, and vegetation. As shown in the comparison (Table 7), the classification based on the UNVI-AT series owns the highest non-crop accuracy. The soil property difference in non-agricultural and agricultural areas can be used to classify the non-crop and crop due to human activities' influence on the soil structure, composition, and trace element contents. For example, cultivation changes to soil organic and elemental contents through chemical  fertilizers, the N b parameter, which represents the background spectral characteristics of soil, is very different between agricultural and non-agricultural areas. The asymmetric double-sigmoid function, the N b is fitted based on the VI value during periods without vegetation cover. The NDVI was calculated using only the RED and NIR bands, and the UNVI was calculated based on all bands, so UNVI was more sensitive to differences among different types of soils. The other advantage of the UNVI is that, unlike the NDVI, its value range is not limited from 0 to 1. Therefore, the UNVI can retain its sensitivity in densely vegetated areas, and differences in its value can reveal differences in vegetation cover. The parameter N a represents the amplitude of change in a VI during the growing cycle. The N a fitted using the UNVI contrasted between noncrop cover, soybean, and corn. But, it was difficult to distinguish non-crop cover using the N a fitted through the NDVI. With a high NDVI value, corn, soybean, and non-crop cover had the same characteristics in N a . The better separability of N a could also improve the separability of the rate of change of the increasing slope, represented by p. Figure 11 shows the separability between agricultural areas and non-agricultural areas based on N a and N b . It was easier to classify agricultural and non-agricultural areas using the UNVI-AT series. As for classification characteristics based on the NDVI-AT series, there is not a distinct boundary between agricultural and non-agricultural example points; it misclassified some soybean fields as the city. It is why the F1 scores of classifying corns based on the NDVI-AT series were lower than that based on the UNVI-AT series.
The advantages of the UNVI also are helpful for the classifier based on the UNVI-AT series to be more accurate at classifying non-crop cover and city areas (Figure 12). The accuracy of classification based on the UNVI-AT series was greater than 90% for non-crop and above 60% for the city. The accuracy of classification based on the NDVI-AT series was only approximately 30% for non-crop and approximately 50% for the city.

Effects of using AT data
We compared the performance of classifications based on the UNVI-AT and UNVI-TIME series. Classification based on the UNVI-AT was more accurate than that based on the UNVI-TIME. Figure 11. Separability between agricultural areas, which include corn and soybean, and non-agricultural areas, which include city areas and non-crop cover, based on the Na and Nb values calculated using the UNVI-AT and NDVI-AT series. Table 8, the difference in classification accuracy, of city, non-crop, soybean, or corn, based on the UNVI-AT series between 2019, 2018, and 2017 was about 5%. However, for classification based on the UNVI-TIME series, the difference in classification accuracy for soybean reached more than 10% . For soybean in 2018, the accuracy of classification was only 71%, compared with that in 2019, the accuracy reduced approximately by 14%. In contrast, there was almost no change in the accuracy of the classification of corn based on the UNVI-AT series.

As shown in
The comparison of the separability between corn and soybean in 2017, 2018, and 2019 between the UNVI-AT series and the UNVI-TIME series is shown in Figure 13. There was relatively little difference between the D d − D i (D d and D i were the parameters of the asymmetric double-sigmoid function) of corn versus soybean from 2017 to 2019, as calculated based on the UNVI-AT series. Although corn could be separated from soybean in 2019 based on the UNVI-TIME series, it was difficult to classify corn and soybean based on D d − D i the UNVI-TIME series in 2018 and 2017. The soybean D d − D i with the UNVI-TIME series in 2018 was larger than that in 2017 and 2019 and very similar to the corn D d − D i .
For non-crop and city cover, there was little change in separability, whether based on the UNVI-AT or UNVI-TIME series (Figure 14). There were apparent differences in the non-crop D d − D i between 2018 and the other two years based on the UNVI-AT series. It was probably led by climate change in 2018. Figure 14 also proved that the UNVI-TIME series might have better performance in classifying city and non-crop because there was litter intra-class difference.  In conclusion, considering temperature characteristics can reduce the differences resulting from changes in the planting period due to environmental factors for corn and soybean classification.

Effects of parameter changes in the asymmetric double-sigmoid function
In most studies, the start and end point of the vegetation-growing period is defined by the extent of change or the derivative change of a VI time series. Although these methods have been effective to  research vegetation phenology characteristics, there were several reasons why these methods have not been used in this research.
First, it is difficult to accurately locate the real start and end points of the vegetation growing period. In most previous studies, the start and end points were defined based on changes in a VI time series. Vegetation phenology can be accurately characterized using VI time series, but not changes in vegetation phenology. It is easier to detect the start and end points using a VI curve.
Second, by locating vegetation phenology change points based on a VI time series, vegetation phenology will change together with the change in the VI time series. Due to the influence of environmental factors, the locations of vegetation phenology points using this method would vary between different years.
In this study, after using the asymmetric double-sigmoid function to fit the relationship between the VI and AT series, we set D s and D e as the start and end points of the vegetation growing period, respectively, by referring to the vegetation phenology change point location method reported by Zhong et al. (2011). The values of D s and D e were calculated as follows: By comparing the D s ( Figure 15A) and D e ( Figure 15B) values calculated using different parameters, it was found that changes of D i , D d , p, and q would make the AT values of D s and D e larger or smaller, but in fact, VI values of D s and D e should not change with the four parameters changing. Although changes of N a did not lead to AT values change of D s or D e , VI values of D s or D e did change when N a changed. Changes in the p and q values could also change the maximum value of the VI. The comparison proved that the higher the value of N a is, the higher the VI values of D s and D e will be. A higher N a value indicates a better vegetation growth state, and the start and end of the growing period, which is calculated based on the VI time series, should not increase under a better vegetation growth state. For the same reason, it is also not logical that VI values of D s and D e calculated using the vegetation phenology function changed with changes in the background N b . At the same time, it was difficult to explain the reason for using the constant term ln 4.5 √ − 2 √ 2 instead of another number because the selection of a start or end point near the located point may also be logical. In our study, we used the difference between the AT with the maximum increase rate of VI, D i , and the AT with the minimum rate of VI decreases, D d , instead of the vegetation growing period D e − D s because the difference of D i and D d was determined by the VI curve and was not influenced by changes in the parameters p or q. Figure 16 shows a comparison of separability between D d − D i and D e − D s based on the UNVI-AT series, where D d − D i has better separability. The worst separability when D e − D s was used was found for corn in 2018 and 2019. The D e − D s of corn in 2018 and 2019 was less than that in 2017. This was probably because the changes in the values of p and q were not only due to the type of vegetation but also due to changes in N a and N b .

The influence factors about the classification of Nebraska
According to classification in Nebraska, the overall accuracy and Kappa of classification changed observably. We calculated the F1-score of each landcover and showed the result in Table 9.
The F1-score of classification proved that the performance of the classifier is stable in terms of corn and soybean. The accuracy of corn and soybean is over 90% all the years. But the grassland and forest accuracy are unstable in three years. The maximum difference of accuracy is over 15%. The difference in accuracy in terms of alfalfa between 2017 and 2018 even exceeded 30%. So, the overall accuracy and Kappa of classification is mainly led by the poor performances of classifications in the forest, grassland and alfalfa.
The poor performance in terms of alfalfa and grassland is led by the unconformity between alfalfa and grassland phenology and Asymmetric Double-sigmoid Function. As shown in Figure 17, there is a big difference between the curves fitted by alfalfa, grassland and city sample points. And the fit curve reduces the difference between alfalfa, grassland and forest.
The alfalfa field is very fragmented and small, so it is very hard to gain sufficient alfalfa sample points. The classifier cannot obtain enough information through alfalfa sample points. So, in the study area, some areas were misclassified as alfalfa. And due to classification trained only with time series, the field of classification was fragmented compared with CDL, and this problem led to some small fields being smaller than that in CDL. And part of fields was mixed with multiple crops.
As for the city, it was mixed with forest and building; it is hard to select pure city sample points by CDL. The city sample point was mixed with building, forest and grassland frequently. So, as shown in Figure 17, the maximum UNVI value of city sample points was about 0.8, which should be 0.2 or lower theoretically. And part of the city was classified as alfalfa, and the accuracy of the city in 2017 and 2018 was lower than 50%.

Conclusion
We built a classifier based on the parameters of the asymmetric double-sigmoid function, which was fitted by a VI time series and AT time series. And the classifier was applied in 2017 without pretraining. We found that the characteristics based on a VI-AT series were stable and unaffected  by changes in environmental factors. According to the crop classification results, the following conclusions could be drawn.
1. The NDVI and UNVI time series could be used to build an asymmetric double-sigmoid relationship with AT for vegetation. The classifier based on the UNVI-AT series had a better classification accuracy than that based on the NDVI-AT series, and the characteristics of the UNVI-AT series were more separable than those based on the NDVI-AT series. 2. Compared with the classifier based on the UNVI-TIME series, the classifier based on the UNVI-AT series maintained a higher classification accuracy from 2017 to 2019. The study proved that the characteristics based on the UNVI-AT series were more stable than those found on the UNVI-TIME series. The time required for vegetation growth changed as environmental factors changed, but the AT required for vegetation growth did not change in the same way.
The study proved that extracting unchangeable vegetation characteristics based on VIs and temperature is practicable, and the accumulated temperature can be a good indicator of vegetation growth. The accuracy of the fitting curve is the foundation of classification accuracy. The asymmetric double-sigmoid function can extract characteristic information in terms of corn, soybean, and forest, but it is unsuitable for extracting characteristics of the city, grassland and alfalfa. So, it is urgent to build up a phenology model, which is suitable for all vegetation landcover in future research.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Data availability statement
The Land Surface Temperature data that support the findings of this study are available in [MOD11A2.006 Terra Land Surface Temperature and Emissivity 8-Day Global 1 km] at [https://developers.google.com/earth-engine/ datasets/catalog/MODIS_006_MOD11A2]; The Sentinal-2 and Landsat-8 image data that support the findings of this study are available at [https://earthexplorer.usgs.gov/]; CDL data that support the findings of this study are available at [https://nassgeodata.gmu.edu/CropScape/].