A Hidden Markov Models Approach for Crop Classification: Linking Crop Phenology to Time Series of Multi-Sensor Remote Sensing Data

Abstract

Vegetation monitoring and mapping based on multi-temporal imagery has recently received much attention due to the plethora of medium-high spatial resolution satellites and the improved classification accuracies attained compared to uni-temporal approaches. Efficient image processing strategies are needed to exploit the phenological information present in temporal image sequences and to limit data redundancy and computational complexity. Within this framework, we implement the theory of Hidden Markov Models in crop classification, based on the time-series analysis of phenological states, inferred by a sequence of remote sensing observations. More specifically, we model the dynamics of vegetation over an agricultural area of Greece, characterized by spatio-temporal heterogeneity and small-sized fields, using RapidEye and Landsat ETM+ imagery. In addition, the classification performance of image sequences with variable spatial and temporal characteristics is evaluated and compared. The classification model considering one RapidEye and four pan-sharpened Landsat ETM+ images was found superior, resulting in a conditional kappa from 0.77 to 0.94 per class and an overall accuracy of 89.7%. The results highlight the potential of the method for operational crop mapping in Euro-Mediterranean areas and provide some hints for optimal image acquisition windows regarding major crop types in Greece.

1. Introduction

The problem of ensuring food security for an increasing population is currently one of the main concerns globally. To solve economic and social issues resulting from current and predicted food shortage, one billion hectares of new cropland would be required in order to meet the demand for food by 2050 [1]. However, taking into consideration environmental restrictions, the potential to expand cropland at the expense of other lands, such as forests or rangelands, is limited [2]. The challenge for agronomists, farmers and their allied partners is to produce humanity’s food in an ecologically sustainable manner, through socially accepted production systems [3]. These trends suggest an increasing demand for dependable and accurate agricultural monitoring to ensure sustainable crop production and investigation of land management practices [4].

Within this framework, spatially explicit cropland information, such as cropland extent and crop type, is crucial to sustain agriculture and preserve natural resources. A basic prerequisite for the implementation of a land-management strategy is the development of up-to-date Land Use/Land Cover (LULC) databases over agricultural landscapes. Indeed, LULC data, regarding the spatial distribution of crop types, is considered as key information from a geostrategic point of view. The research community is moving towards providing and timely agricultural maps at national or global level of detail [5].

Over the past decade satellite images offer a valuable source of information concerning the monitoring of the Earth’s surface in fine spectral and spatial scales. Satellite based earth observation has been used to map crop types under a variety of environmental conditions, providing synoptic coverage of fields in several spectral regions, smooth integration with existing geographical databases under a cost-effective and time-saving approach than traditional statistical surveys [6]. Through these studies, remote sensing techniques have proven to be cost-effective in widespread agricultural lands in Africa, America, Europe and Australia.

Monitoring and mapping vegetation involves investigating vegetation dynamics, such as phenological states and the seasonal growth of crop types. Spectral behavior of agricultural parcels is constantly changing; different crop types may be at a certain instant in the same phenological state, depicting similar spectral attributes, but diverge remarkably in another instant. Classification of parcels based on single-date images, even if they are acquired at critical growth states, cannot offer reliable results in the case of crops with similar growing cycle. As a result, the significance of classification based on multi-temporal images has been well-recognized and especially regarding vegetation mapping, the usage of seasonal imagery is vital [5,7,8,9,10,11,12].

Low resolution images with high revisit frequency have been processed on the continental and global scale, providing consistent information at high temporal resolution while covering large areas at low costs [13]. However, because of sub-pixel heterogeneity, the spatial resolution of the imagery may result in significant errors in the estimated crop areas [14,15]

At the regional level, crop area estimations have been significantly improved since the introduction of the MODIS sensor with 250 meter ground resolution [13]. MODIS offered unprecedented capabilities for large-area LULC mapping by providing global coverage, half-day revisit capacity and intermediate spatial resolution. Several studies have already demonstrated successfully the potential of these data for detailed LULC mapping in an agricultural setting [15], especially for areas where the typical field size is large [5,16].

Multi-spectral, medium-high resolution images, acquired mainly by Landsat Thematic Mapper/Enhanced Thematic Mapper Plus (TM/ETM+), SPOT and RapidEye [17,18,19,20], have been used for regional to local scale crop mapping, either on mono-temporal basis or under a multi-temporal perspective, accounting for small sized fields or heterogeneous crop patterns. To solve the issue of high within-class variability originating from various agronomic practices, the synergy of multi-temporal optical and Synthetic Aperture Radar (SAR) data has also been proposed; the increased set of multi-temporal imagery enables the continuous monitoring of all stages of vegetation development [21]. While multi-temporal data of medium-to-high resolution offers high potential for crop discrimination in fine-structured agricultural landscapes, integration of the temporal information in the classification process is not trivial. Precise annual mapping of crops, through an approach that could be used routinely over large areas, remains challenging [22,23].

Although a variety of algorithms has been employed in crop mapping studies, including among others, the minimum distance, Mahalanobis distance, maximum likelihood, spectral angle mapper and support vector machines, an approach integrating phenological models in the classification process has been given little attention so far. Through phenology, remote sensing observations and biophysical changes during vegetation’s growth can be linked statistically in order to discriminate crop types. A possible way of incorporating knowledge of phenology into the classification process lies on the adoption of the stochastic Hidden Markov Models (HMMs). HMMs allow the simulation of crop dynamics, exploiting the spectral information of their phenological states and their relations. In this regard, a common assumption is made that the vegetation signal and the different phenological states are considered random variables [24]. The correlation of phenological states is described by different transition state probabilities in each crop model. As far as the different cultivation practices, the algorithm reckons the possibility of temporal variation in the phenological cycle. Different growing states of the same crop type per image can be introduced instead of using a generalized crop model. These states correspond to different spectral attributes and are used jointly to define each model. This is the basic advantage of HMMs compared to other techniques that produce simulations of average seasonal phenology [25] and may fail to account for a restrained or accelerated phenological progress.

Previous work has tested the application of HMMs in Landsat time series to classify mountain vegetation in Norway [26] and arable land in Brazil [27], in MODIS-NDVI time series covering cultivated areas of the United States [28] and NDVI data derived from the Advanced Very High Resolution Radiometer (AVHRR) over the West African savanna [24]. In all the aforementioned studies, the low and medium resolution images have been reported to be adequate for the classification of large-sized agricultural holdings. However, Mediterranean regions that are characterized by distinct environmental and climate settings, high spatio-temporal ecological heterogeneity [29,30], variety of crop types and high fragmentation of farming lands [23,31], require a different approach.

The main aim of this study is the development of a robust crop mapping technique adopting the theory of Markov chains and phenological models, over a Euro-Mediterranean agricultural area. Specifically, this work proposes a crop classification approach that integrates high and medium resolution remote sensing images in order to monitor constant variations in the ecological process of the cropping systems. A pixel-based methodology was selected, instead of using the segment-based approach proposed by [27], to avoid errors produced by segmentation algorithms. The per-segment approach applied to small-sized crop parcels, found over Euro-Mediterranean areas, may have the disadvantage of falsely including within field-crop objects small non-vegetated classes (i.e., roads, canals) leading to an overestimation of the total vegetated area [32].

In particular, the objectives of the study are: (1) the identification of different crop types using a sequence of four seasonal multispectral Landsat ETM+ and a RapidEye image, processed simultaneously through Hidden Markov Models, (2) the assessment of the impact on the accuracy of a pan-sharpening procedure applied to the lower resolution Landsat ETM+ imagery and (3) the investigation of the role of the temporal resolution and extent of the image sequence used, in relation to the phenological cycle of each crop type. The multi-sensor and multi- temporal approach is motivated by the acknowledgment of the potential of coarser spatial resolution data to cover large geographic extends, the demands of complex territories and the growing interest in exploiting multi-scale data synergistically [4]. Furthermore, the definition of optimal temporal acquisition windows is considered vital by several crop mapping studies [8,11,18,33] while it can improve classification accuracy significantly.

2. Study Area

The research site is an irrigated agricultural area, near the city of Thessaloniki, Greece. The study area is dominated by rice and cotton while maize, sugar beet, wheat and alfalfa are planted to a smaller extent. The cropping calendar (planting and harvesting dates) of the area’s crop types is presented in Figure 1.

Figure 1. Idealized cropping calendar of the main crop types grown in the study area.

Figure 1. Idealized cropping calendar of the main crop types grown in the study area.

Rice, cotton, maize and sugar beet are summer crops and those fields are characterized by dense vegetation during summer months. Wheat is harvested before June and is a spring crop. Rice, cotton, maize, sugar beet and wheat are considered annual crops and have a 12-months cycle. Alfalfa on the other hand can have 3–4 cuttings and flowerings per year, usually between May and September. Thus, the cropping pattern of the study area can be considered heterogeneous regarding the dates of planting, emergence, and harvesting. The majority of the parcels are rectangular but small sized. Despite the applied land consolidation the size of the parcels ranges from 0.006 to 10 ha. The terrain across the study area is relatively flat. The average annual temperature of the study area is 15.8 °C. The area is characterized by modest annual rainfall, averaging 441 mm/y.

3. Materials and Methods

3.1. Outline of the Methodology

Originally, the Landsat ETM+ images were registered to the higher resolution RapidEye image, which has been georeferenced using ground control points (GCPs) identified over VHR orthophotographs (Figure 2) in the same geodetic system with the vector dataset representing field entities of the area (Land Parcel Identification System-LPIS). LPIS was visually corrected for small inconsistencies. The multispectral ETM+ images were pan-sharpened using the panchromatic band and the High Pass Filter (HPF) algorithm. Four synthetic images were produced with a spatial resolution of 15 meters. Nine different classifications experiments were applied on image sequences with variable spatial and temporal characteristics. A common set of training data, derived from the LPIS, was used to estimate the parameters of the crop models. For each HMM, we calculated the probability that the specific set of temporal observations corresponded to a class. Each pixel is assigned to the class whose crop-model emits the maximum probability. The results of the classification tests were evaluated in terms of overall accuracy and kappa coefficients.

Figure 2. Overall process diagram of this research study.

Figure 2. Overall process diagram of this research study.

3.2. Satellite Data and Preprocessing

Imagery acquired from sensors with different spectral, spatial and radiometric characteristics was used in the analysis. More specifically, the multispectral (excluding the thermal bands) and panchromatic components of four Landsat ETM+ images (184/32) and one multispectral RapidEye image, all acquired on different dates of 2010 (Figure 3 and Figure 4), were employed in this study (Table 1).

Table 1. Description of the satellite data used in the study. *: blue; G: green; R: red; NIR: near-infrared; SWIR: shortwave-infrared bands.

Time Step	Sensor	Date of Acquisition	Spatial Resolution Multi/Pan	Radiometric Resolution	Spectral Bands *
t1	Landsat ETM+	07/05/2010	30/15 m	8-bit	B, G, R, NIR, SWIR1, SWIR2
t2	Landsat ETM+	08/06/2010	30/15 m	8-bit	B, G, R, NIR, SWIR1, SWIR2
t3	RapidEye	05/08/2010	5 m	16-bit	R, G, B, Red edge, NIR
t4	Landsat ETM+	27/08/2010	30/15 m	8-bit	B, G, R, NIR, SWIR1, SWIR2
t5	Landsat ETM+	30/10/2010	30/15 m	8-bit	B, G, R, NIR, SWIR1, SWIR2

Table 1. Description of the satellite data used in the study. *: blue; G: green; R: red; NIR: near-infrared; SWIR: shortwave-infrared bands.

Time Step	Sensor	Date of Acquisition	Spatial Resolution Multi/Pan	Radiometric Resolution	Spectral Bands *
t1	Landsat ETM+	07/05/2010	30/15 m	8-bit	B, G, R, NIR, SWIR1, SWIR2
t2	Landsat ETM+	08/06/2010	30/15 m	8-bit	B, G, R, NIR, SWIR1, SWIR2
t3	RapidEye	05/08/2010	5 m	16-bit	R, G, B, Red edge, NIR
t4	Landsat ETM+	27/08/2010	30/15 m	8-bit	B, G, R, NIR, SWIR1, SWIR2
t5	Landsat ETM+	30/10/2010	30/15 m	8-bit	B, G, R, NIR, SWIR1, SWIR2

Landsat ETM + sensor launched in April 1999 has a spatial resolution of 30 meters for the six reflective bands, 60 meters for the thermal band, and 15 meters for the panchromatic (pan) band. On 31 May 2003, the ETM+ Scan Line Corrector (SLC) failed causing the scanning pattern to exhibit wedge‐shaped, scan-to-scan gaps, which are most pronounced along the edge of the scene. The scans give near-contiguous coverage of the surface scanned below the satellite in the center of the image (approximately 22 km wide).

RapidEye is a constellation of 5 multispectral satellite sensors launched in August 2008 with a primary focus on agricultural applications. The RapidEye sensor has a multispectral push broom imager with a spatial resolution of 6.25 meters. It captures data in five spectral bands covering visible–infrared part of the electromagnetic spectrum: blue (440–550 nm), green (520–590 nm), red (630–685 nm), red edge (690–730 nm), and near infrared (760–850 nm). In our study, we used the RapidEye (Level 3A) in which radiometric, sensor, and geometric correction have been applied and resampled to a 5 meters spatial resolution. We georeferenced the RapidEye imagery to the Greek Geodetic Reference System 1987, using ground control points, identified on natural color orthoimages with 50 cm spatial resolution acquired on 2007. Between 20 and 30 ground control points were used to co-register the Landsat scenes to the RapidEye imagery, with a Root Mean Square error (RMS) of less than a pixel. We did not apply any atmospheric correction or radiometric normalization since the adopted classification approach does not employ any direct comparison of pixel DN values between the temporal sequences of images. Instead, the classification scheme assigns pixels to crop classes, according to their similarity of states within each image separately and in a subsequent step the temporal images are linked using statistical relationships. In this context, atmospheric correction was not considered a prerequisite [32].

The same geographic subset was identified on every Landsat imagery, with no effective clouds or sensor defects, such as the “SLC-off problem”, covering approximately 7500 ha of cultivated area (1760 by 1710 RapidEye pixels). Additionally, regarding the Landsat images, the panchromatic images were merged with the multispectral ones, using the High Pass Filter (HPF) algorithm and four synthetic images were produced with a spatial resolution of 15 meters [21]. Finally, in order to achieve spatial correspondence for each pixel, all Landsat ETM+ images (original and pan-sharpened) were re-sampled to 5 meters, using a nearest-neighbor algorithm to match the spatial resolution the RapidEye image (Figure 3).

Figure 3. The set of images used included four Landsat ETM+ images (image t₁, image t₂, image t₄, image t₅) illustrated in false-color composite (R: NIR, G:Red, B:Green) and one Rapideye image (image t₃) in false-color composite (R: NIR, G:Red, B:Green).

Figure 3. The set of images used included four Landsat ETM+ images (image t₁, image t₂, image t₄, image t₅) illustrated in false-color composite (R: NIR, G:Red, B:Green) and one Rapideye image (image t₃) in false-color composite (R: NIR, G:Red, B:Green).

Figure 4. The acquisition dates of the images used in the study.

Figure 4. The acquisition dates of the images used in the study.

3.3. Reference Data

The Land Parcel Identification System (LPIS) is a fundamental part of the Integrated Administration and Control system that has been developed and adopted in 1992 by the EU as the spatial component for the implementation and supporting of the Common Agricultural Policy (CAP) and land management across Europe [34]. The main functions of the LPIS are localization, identification and quantification of agricultural land via very detailed geospatial data, in order to spatially represent the activities of farmers on their land and facilitate the geographical identification of the agricultural parcels declared annually to receive funding [35].

Although the regulatory requirements for the LPIS are uniform across the EU sector, the particular implementations are subject to member states. The Greek GIS-based LPIS integrates information about the crop type, the acreage of a parcel, the identity of the farmer and relates it to a vector layer comprised of the declared parcels. Since the information of this database is gathered from the declaration of the farmers it cannot be considered flawless. In this respect, an expert from “Greek Payment Authority of Common Agricultural Policy Aid Schemes of the Ministry of Rural Development and Food” visually examined and corrected the parcels’ crop type and boundaries manually, taking into account the cropping calendar of the study area and the spectral- temporal profile of each parcel. The detailed delineation of the boundaries was guided mainly by the high resolution RapidEye image. In total, 3319 declared parcels were found in the study area. A set of 55 parcels was used as training data and the rest was used during the accuracy assessment of the classification.

3.4. Description of the Proposed HMMs Classification Algorithm

The temporal evolution of vegetation can be described effectively by the state-oriented approach of Hidden Markov Models (HMMs). Each cultivated parcel has a dynamic behavior that depends on cropping phenology, climatic conditions, drought, water irrigation and chemical nutrients (Figure 5).

Figure 5. Temporal sequence of images t₁, t₂ and t₃ covering the same area containing parcels cultivated by various crop types (1 = maize, 2 = rice, 3 = wheat, 4 = alfalfa and 5 = sugarcane). It is indicated that the different phenological states of crops at each time step impose significant variance in the between-class separability.

Figure 5. Temporal sequence of images t₁, t₂ and t₃ covering the same area containing parcels cultivated by various crop types (1 = maize, 2 = rice, 3 = wheat, 4 = alfalfa and 5 = sugarcane). It is indicated that the different phenological states of crops at each time step impose significant variance in the between-class separability.

Figure 6. Different subsets of satellite images t₃ and t₄ containing parcels with various crop types (1 = cotton, 2 = alfalfa and 3 = maize). In the left subset it can be observed that certain parcels of cotton and alfalfa may resemble according to their phenological state, while other parcels of the same crops can have distinct spectral properties. This can be also observed in the right subset referring to maize fields with different spectral characteristics.

Figure 6. Different subsets of satellite images t₃ and t₄ containing parcels with various crop types (1 = cotton, 2 = alfalfa and 3 = maize). In the left subset it can be observed that certain parcels of cotton and alfalfa may resemble according to their phenological state, while other parcels of the same crops can have distinct spectral properties. This can be also observed in the right subset referring to maize fields with different spectral characteristics.

Furthermore, neighboring parcels of the same crop type, over the same time step, may be at different states of growth due to varying agronomic practices; different planting or harvesting dates and fertilizers can accelerate or restrain the phenological progress (Figure 6).

Given that each parcel changes constantly from state to state (Figure 5) and that each state cannot be directly linked to a remote sensing measurement but to a probability distribution of observations, an HMM can be used to simulate the cycle of vegetation based on statistical relations. In this case, an HMM is a doubly embedded stochastic process comprised by two chains: the external chain of the remote sensing observations and the internal chain of states, which are unknown [24] (Figure 7).

Figure 7. Schematic description of the basic elements of the proposed HMMs.

Figure 7. Schematic description of the basic elements of the proposed HMMs.

In this study, a crop classification algorithm, based on Markov-chain analysis, was designed and implemented using Matlab^®. Accordingly, a model was built describing the phenological states during the dates of the study for each crop type. The sequence of observations consisted of a set of remote sensing measurements O = {O_t1 ,..., O_tn}, where t = {t₁ ,..., t_n} are the acquisition dates of the images. The hidden states correspond to the different phenological states S = {S₁ ,..., S_m} of each crop type and Q = {q_t1,..., q_tn} is the fixed sequence of hidden states. During the acquisition of the satellite data, the identified states of the parcels in the study area were: S₁, no vegetation or beginning of emergence, S₂, medium vegetation or growth state, S₃, dense vegetation or flowering and S₄, dry vegetation or harvesting state.

An HMM is characterized by the following elements:

• The state transition probability matrix A, where a_i,j(t) = P[q_t = S_j |q_t-1 = S_i], denotes the transition probability from state S_i to state S_j at time t. The emission probability matrix B, defines the probability that O_t is emitted by state S_j, i.e., b_j(O_t) = P[O_t | q_t = S_j ]. In order to estimate the symbol probability distributions B, a multivariate Gaussian distribution is assumed for the observed spectral data. The mean vector ${\mu }_i$ and the covariance matrix ${\sum }_i$ were calculated by the training data for each crop type, for each state and for every image by the equation,

  $$b_i\left(O_t\right)=\frac{1}{\sqrt{{\left(2\pi \right)}^{d }|{\Sigma }_i|}}\exp \left(-\frac{{\left(O_t-{\mu }_i\right)}^T{\Sigma }_i^{-1}\left(O_t-{\mu }_i\right)}{2}\right)$$

  (1)
• The initial probability π_i is the probability of being in state Si at time t1, i.e., π_i = P[q_t1 = S_i]. The parameters A, B and π_i were estimated by the set of training data. The set of training samples was selected and defined according to our knowledge of local agronomic practices and the cropping calendar of the district. To ensure classification success, all classes need to be described by representative training samples. The samples define the different states of crop types according to their different spectral attributes. The set of the training data was evenly distributed in the study area, located in homogenous fields and not in boundary mixed pixels. It should also be noted that in each image, fields of the same crop type may be in different states; usually a state before or a state after (Figure 6). Judging by the cropping calendar and the dates of the used images, not all transitions between states are possible in this case study. Once we have estimated the parameters A, B, π_i of each model λ and given the sequence of observations O for each pixel, the probability that the sequence O was generated by these models λ is defined by:

  $$\text{P}\left({\text{q}}_1={\text{S}}_{\text{i}},\dots ,{\text{q}}_{\text{t}}={\text{S}}_{\text{j}},{\text{O}}_1,\dots ,{\text{O}}_{\text{t}}\text{|λ}\right)=\text{P}\left({\text{O}}_1,\dots ,{\text{O}}_{\text{t}}{\text{ |q}}_1={\text{S}}_{\text{i}},\dots ,{\text{q}}_{\text{t}}={\text{S}}_{\text{j}}\right)\times \text{P}\left({\text{q}}_1={\text{S}}_{\text{i}},\dots ,{\text{q}}_{\text{t}}={\text{S}}_{\text{j}}\text{ |λ}\right)={\text{π}}_{{\text{q}}_1}\times {\displaystyle \prod }_{\text{l}=2}^{\text{t}}{\text{a}}_{{\text{q}}_{\text{t}-1,}{\text{q}}_{\text{t}}}\left(\text{l}\right)\times {\displaystyle \prod }_{\text{l}=1}^{\text{t}}{\text{b}}_{{\text{q}}_{\text{t}}}\left({\text{O}}_{\text{l}}\right)$$

  (2)

Detailed mathematical explanation of HMMs has been reported in previous studies [24,36] which propose the implementation of the Forward algorithm to simplify computations of Equation (2).

In this study, five models were built (one for each crop class), and five measurements of probability were estimated for each pixel (Equation (2)). Let us consider a temporal spectral sequence of pixel x belonging to l_m, the most probable crop, where L = {l₁,..., l_k} is the set of possible crop classes. For each pixel x, the most likely crop-class is determined by the following rule:

$$x\in {\lambda }_m when {\lambda }_m=\arg \max \left(\text{P}\left({\text{q}}_{\text{1}}{=\text{S}}_{\text{i}}{\text{,…, q}}_{\text{t}}{=\text{S}}_{\text{j}}{\text{, O}}_{\text{1}}{\text{,…, O}}_{\text{t}}{\text{|λ}}_{\text{m}}\right)\right)$$

(3)

Finally, nine different HMM models were developed in order to evaluate the role of the temporal resolution and extent of the image sequence used, in relation to the phenological cycle of each crop type, as well as to assess the utility of the pan-sharpening procedure in terms of classification accuracy (Table 2).

Table 2. Selections of different images included in the evaluation of our proposed methodology.

	Imagery Employed									Rationale of the Classification Experiment
	t1	t2	t4	t5	t1	t2	t4	t5	t3
	Original ETM+				Pan-Sharpened ETM+				RapidEye
HMM-1	Χ	Χ	Χ	Χ					Χ	Assessment of the HMMs approach
HMM-2					Χ	Χ	Χ	Χ	Χ	Influence of pan-sharpening
HMM-3					Χ	Χ	Χ	Χ		Influence of spatial resolution-RapidEye
HMM-4					Χ	Χ			Χ	Influence of temporal extent and resolution
HMM-5						Χ		Χ	Χ	Influence of temporal extent and resolution
HMM-6						Χ	Χ		Χ	Influence of temporal extent and resolution
HMM-7					Χ				Χ	Influence of temporal extent and resolution
HMM-8					Χ		Χ			Influence of temporal extent and resolution
HMM-9						Χ	Χ			Influence of temporal extent and resolution

Table 2. Selections of different images included in the evaluation of our proposed methodology.

	Imagery Employed									Rationale of the Classification Experiment
	t1	t2	t4	t5	t1	t2	t4	t5	t3
	Original ETM+				Pan-Sharpened ETM+				RapidEye
HMM-1	Χ	Χ	Χ	Χ					Χ	Assessment of the HMMs approach
HMM-2					Χ	Χ	Χ	Χ	Χ	Influence of pan-sharpening
HMM-3					Χ	Χ	Χ	Χ		Influence of spatial resolution-RapidEye
HMM-4					Χ	Χ			Χ	Influence of temporal extent and resolution
HMM-5						Χ		Χ	Χ	Influence of temporal extent and resolution
HMM-6						Χ	Χ		Χ	Influence of temporal extent and resolution
HMM-7					Χ				Χ	Influence of temporal extent and resolution
HMM-8					Χ		Χ			Influence of temporal extent and resolution
HMM-9						Χ	Χ			Influence of temporal extent and resolution

3.5. Accuracy Assessment

Each classification test was evaluated in terms of overall accuracy (OA) and Kappa coefficient, by comparing the reference data with the classified images, pixel by pixel. Overall accuracy represents the proportion of the correctly classified pixels relative to the total number of validation pixels. Kappa coefficient takes into account all the elements of the error matrix and is a measure of the proportional improvement by the classifier over a purely random assignment to classes [37]. Despite the fact that both the OA and Kappa coefficient measure the agreement between the classified map and the reference data, Kappa is often considered a better indicator of classification performance because it excludes chance agreement [38]. Finally, the conditional Kappa coefficient was used for assessing the agreement for the individual crop categories of the maps.

4. Results and Discussion

The results of the accuracy assessment of the classification tests are presented in Table 3 including overall accuracy, overall kappa coefficient and kappa coefficient of each class and are further discussed in the following sections.

Table 3. Overall accuracy, overall kappa coefficient and the kappa coefficient of each class for all experiments according to the selected images.

		HMM-1	HMM-2	HMM-3	HMM-4	HMM-5	HMM-6	HMM-7	HMM-8	HMM-9
Overall Accuracy		84.7%	89.7%	88.5%	87.5%	90.5%	92.5%	75.9%	76.0%	91.1%
Overall Kappa coefficient		0.774	0.843	0.825	0.811	0.852	0.881	0.658	0.655	0.859
Conditional Kappa Coefficient	Cotton	0.662	0.770	0.743	0.722	0.823	0.894	0.504	0.570	0.871
	Rice	0.911	0.882	0.875	0.848	0.924	0.907	0.730	0.640	0.898
	Sugar beet	0.954	0.936	0.934	0.930	0.922	0.906	0.875	0.854	0.900
	Alfalfa	0.789	0.852	0.897	0.795	0.606	0.788	0.694	0.858	0.883
	Maize	0.742	0.861	0.831	0.859	0.867	0.805	0.834	0.641	0.702
	Wheat	0.930	0.916	0.869	0.953	0.659	0.690	0.895	0.705	0.681

Table 3. Overall accuracy, overall kappa coefficient and the kappa coefficient of each class for all experiments according to the selected images.

		HMM-1	HMM-2	HMM-3	HMM-4	HMM-5	HMM-6	HMM-7	HMM-8	HMM-9
Overall Accuracy		84.7%	89.7%	88.5%	87.5%	90.5%	92.5%	75.9%	76.0%	91.1%
Overall Kappa coefficient		0.774	0.843	0.825	0.811	0.852	0.881	0.658	0.655	0.859
Conditional Kappa Coefficient	Cotton	0.662	0.770	0.743	0.722	0.823	0.894	0.504	0.570	0.871
	Rice	0.911	0.882	0.875	0.848	0.924	0.907	0.730	0.640	0.898
	Sugar beet	0.954	0.936	0.934	0.930	0.922	0.906	0.875	0.854	0.900
	Alfalfa	0.789	0.852	0.897	0.795	0.606	0.788	0.694	0.858	0.883
	Maize	0.742	0.861	0.831	0.859	0.867	0.805	0.834	0.641	0.702
	Wheat	0.930	0.916	0.869	0.953	0.659	0.690	0.895	0.705	0.681

4.1. Multitemporal Classification Using the Original ETM+ and RapidEye Imagery

It has been observed [6] that the spatial resolution of the imagery should be at or below the size of the fields. Nevertheless, detailed information provided by high resolution images does not meet the requirements for temporal availability and cost-effective processing framework. When lower resolution sensors are selected, the accuracy of the classification is affected by the mixed pixel problem. For this reason, we explored the processing of time series of high and medium resolution images simultaneously. In the first classification experiment (HMM-1), considering one RapidEye and the original multispectral Landsat ETM+ images (Table 2), the classification model achieved an overall accuracy of 84.7% and an overall Kappa coefficient of 0.774. The conditional kappa coefficients ranging from 0.662 to 0.954, suggest that the spectral information was adequate for discriminating the majority of the crop types.

4.2. Multitemporal Classification Using the Pan-Sharpened ETM+ and RapidEye Imagery

In order to evaluate the contribution of the pan-sharpening of the Landsat ETM+ images in the improvement of the classification process, the developed HMM was tested using the synthetic ETM+ bands along with the multispectral RapidEye imagery (HMM-2). As the spatial resolution of the Landsat ETM+ images increased, the performance of the classification model improved, reaching 89.7% in overall accuracy and 0.843 in overall kappa coefficient (which corresponds to an increase of 5% and 0.069 respectively).

Figure 8. Visual assessment of the classification results obtained from experiment HMM-1 considering the original ETM+ and RapidEye images (a) and from experiment HMM-2 considering the pan-sharpened ETM+ and RapidEye images (b). HMM-1 experiment resulted to more extended classification errors along the parcel’s borderline (c) and the respective classification errors (d) compared to lower classification confidence (e) and confidence score (f) obtained from experiment HMM-2.

Figure 8. Visual assessment of the classification results obtained from experiment HMM-1 considering the original ETM+ and RapidEye images (a) and from experiment HMM-2 considering the pan-sharpened ETM+ and RapidEye images (b). HMM-1 experiment resulted to more extended classification errors along the parcel’s borderline (c) and the respective classification errors (d) compared to lower classification confidence (e) and confidence score (f) obtained from experiment HMM-2.

Regarding the conditional kappa coefficients of individual crops, the highest increase of ~0.110 is observed for “cotton” and “maize” classes, which presented the lowest discrimination ability in the previous classification experiment (HMM-1). This sharp increase can be attributed to the shape of the respective crop fields, being more elongated and narrow compared to the other crop fields of the area. The spatial explicit assessment of the classification errors distribution resulting from the HMM-1 and HMM-2 maps (Figure 8), verifies the contribution of the pan-sharpening procedure in the improvement of the classification result. Differences are observed along the boundaries of the different crop fields with a larger proportion of misclassified area found along the borderline of the fields in the case of the original dataset (HMM-1). In addition, the confidence score derived by the computed probabilities of the HMMs and the corresponding confidence images (Figure 8) prove that a larger proportion of each agriculture field is classified with higher confidence in the case of the HMM-2 model.

4.3. Multitemporal Classification Using the Pan-Sharpened ETM+ Imagery

To quantify the contribution of the RapidEye image in the classification scheme the third classification experiment (HMM-3) included only the pan-sharpened Landsat ETM+ images. The overall accuracy of this experiment was 88.5% while the overall kappa coefficient was 0.825. Compared to the results obtained from the previous experiment, integrating the RapidEye imagery, a decrease in both accuracy metrics is evident (1.2% in OA and 0.018 in Kappa values); this can be attributed to the information content inherent to the higher spatial resolution RapidEye imagery.

Parcels not being wide enough to be mapped in a 15 meter resolution may be distinguished in a RapidEye image of 5 meters pixel size. Thus, the RapidEye image actually adds information involving narrow parcels and boundary pixels that cannot be viewed in a Landsat ETM+ scene. Comparison of the individual class results obtained by classification experiments HMM-2 and HMM-3 indicates that the lowest differences exist for the conditional kappa coefficient of class “sugar beet”. This relates to the fact that these crop fields within the study area, have a mean size of 1.28 ha. Their relatively large size allows satisfactory discrimination despite the coarse spatial resolution of the pan-sharpened ETM+ bands.

4.4. Multitemporal Classification Considering Different Temporal Extents

One of our objectives was to assess classification accuracy obtained by decreasing the number of images used in our classification model. Usually five images per year are used to perform multi- temporal classification in agriculture applications [12]. During classification experiments of HMM-4 to HMM-9, the number of images employed in the model decreases, but the resulting overall accuracy does not decrease proportionally. The selection of different temporal images affects the discrimination of certain crop types. By using three images (HMM-4, HMM-5 and HMM-6) instead of five, the conditional kappa coefficient of each summer crop remains relatively high (above 0.722). As a matter of fact, in HMM-6, where the three selected images were all acquired during summer, the corresponding summer crops attain the highest accuracy values ranging between 0.805 and 0.907. When using just three or two images neither the overall accuracy nor the overall kappa coefficient necessarily decrease, but certain individual crop accuracies fall below average. In the cases of the classification experiments HMM-5, HMM-6 and HMM-9, which do not incorporate the ETM+ image acquired on May (t₁), the overall accuracy increases, ranging between 90.5% and 92.5% with an overall kappa coefficient between 0.852 and 0.881. However, these classification experiments perform poorly for class “wheat” because it is an “early crop”. This implies that only in the early May image (t₁) “wheat” fields appear vegetated. In fact, the kappa coefficient of “wheat” drops at 0.659–0.690, while it reaches values of 0.869–0.953 in all other tests. As for class “alfalfa”, it has been already stated that it has a different phenological cycle, depending on the varied dates of the cuttings. This crop type does not have a seasonal pattern like the other crops because the state of “emergence” can be repeated 3–4 times per year and even neighboring fields can differ spectrally. For this crop type at least 4 images are required to reach a stable classification result, as indicated by the classification results.

Comparison of all the classification experiments, suggests that HMM-2 involving four pan-sharpened ETM+ and one RapidEye images, provided the best results, as far as crop-specific accuracies are concerned. It ranges between 0.770 and 0.936, whereas in other tests it falls below 0.60. The visual assessment of the classification errors and the confidence images, verifies the findings of the accuracy assessment, and highlights the importance of the spatial information inserted into this model through the pan-sharpening procedure for improving the classification of the borderline pixels. The pan-sharpening procedure also improved classification of pixels located within narrow fields by better preserving their spectral attributes. Even though, the computed measures of 89.7% in the overall accuracy and 0.843 in the overall kappa coefficient are not the highest, all accuracies per class reach satisfactory levels. Similar performance metrics were reached in [27] and [18], where five images were sufficient to reach an overall accuracy of 85%–89%. The further decrease in the number of images had a significant impact deteriorating the average accuracy per class. As far as the kappa coefficients are concerned, the highest values are observed on class “sugar beet” and class “wheat” (0.936 and 0.916 respectively). This is due to the size of the sugar beet parcels which can be monitored by the Landsat ETM+ scenes. The high accuracy of “wheat” is justified by its different phenology compared to the rest of the crop types. The poorest results were obtained for class “cotton” (0.77 in the kappa index), which was confused with the other “summer crops” due to its high internal spectral variability, resulting from the different farming practices applied. Apart from errors related to sub-pixel heterogeneity, classification errors related to whole field misallocation might arise (Figure 8). This kind of errors stem from various agronomic practices (i.e., different dates of planting and harvesting, usage of herbicides and fertilizers, etc.). This is a common problem in crop mapping studies [21] and could be resolved with the insertion of additional images representing more adequately high within-class variability.

In conclusion, the results suggest optimal dates of scenes according to which crop types are to be examined. It should be noted that deciding on the optimal number and dates of scenes depends on the study area, the number of classes and the variety of cropping systems. In [11] it was proposed that four scenes were adequate for the classification of six crop classes in Japan, while in [32] 2–3 multi-date images were used to discriminate seven classes including “grassland” in the Netherlands. In this paper, by selecting 3–5 scenes the overall accuracy ranges between 87.5% and 92.5%, the overall kappa coefficient between 0.811 and 0.881 and the least conditional kappa coefficient is about 0.60. However, reducing the number of images to two, led to an overall accuracy and overall kappa coefficient of about 75% and 0.65 respectively and a conditional kappa index of 0.50 which are considered moderate. Yet, if we need to constrain to using only two multi-date images we should choose a combination of May and August to account for both “winter” and “summer” cropping systems.

5. Conclusions

In this study, the classification approach was directed to meet the needs of a Mediterranean agricultural area. The approach integrated the following ideas: (a) the theory of HMMs to describe the dynamics of vegetation, (b) the combination of multi-sensor data and (c) the implementation of image enhancement techniques.

The challenge of crop mapping stems from the variety of cropping systems, distinct climate settings and cultivation practices. By using Hidden Markov Models we were able to set a dynamic model per crop type to represent the biophysical processes of agricultural land. Due to the high fragmentation of the land a multi-resolution approach was introduced. Experimental results demonstrated that the integration of even one very high resolution image and pan-sharpening of the set of Landsat images improved the overall classification accuracy by 1.2% and 5% respectively. Spatial explicit classification errors demonstrated that our methodology succeeded in mapping even small-sized fields enhancing classification even on the boundaries of different crop fields. It is worth mentioning that our approach can incorporate and process different satellite data without the implementation of atmospheric correction because the covariance matrix and the mean vector are computed by the samples for each image separately.

However, an evident shortcoming is that models cannot be directly transferred to another year or a different region. This can be attributed to two factors: the inter-annual variation of climate conditions and the heterogeneity of cropping patterns of distinct territories. Even within similar agricultural areas, crops may grow in a different rate depending on the soil or irrigation system. To this end, the algorithm can be extended to integrate weather information and improve the estimation of the transition probabilities. Another drawback is that ground reference data are required each year and over the specific study area to train the classifier. In this paper, to avoid time consuming and costly on-the-field visits, we proposed the use of the ancillary crop maps.

Depending on the application and the investigated crop types, different sets of temporal images may be used. If we are interested in mapping only “summer” or “winter” crops a set of only two images may provide adequate classification results. However, it has been shown that for more complex agricultural landscapes a robust classification scheme requires at least 4–5 images to achieve a kappa index per class above 0.70.

The increased availability of imagery by recent and up-coming satellite missions (i.e., Landsat-8 and Sentinel-2) will offer a more dense set of observations and will broaden the mapping capabilities of the proposed technique. Additionally, further research may explore the integration of vegetation indices (VIs) to examine the impact on the classification performance. Considering the computational complexity of the algorithm, a reduction in the dimension of the input data by using VIs will improve the efficiency of the processing. In the future, classification errors could be eliminated by extending the model to incorporate contextual knowledge and accounting for possible interactions of neighboring pixels.