Particle Filter Approach for Real-Time Estimation of Crop Phenological States Using Time Series of NDVI Images

Abstract

Knowing the current phenological state of an agricultural crop is a powerful tool for precision farming applications. In the past, it has been estimated with remote sensing data by exploiting time series of Normalised Difference Vegetation Index (NDVI), but always at the end of the campaign and only providing results for some key states. In this work, a new dynamical framework is proposed to provide real-time estimates in a continuous range of states, for which NDVI images are combined with a prediction model in an optimal way using a particle filter. The methodology is tested over a set of 8 to 13 rice parcels during 2008–2013, achieving a high determination factor $R^2=0.93$ ($n=379$) for the complete phenological range. This method is also used to predict the end of season date, obtaining a high accuracy with an anticipation of around 40–60 days. Among the key advantages of this approach, phenology is estimated each time a new observation is available, hence enabling the potential detection of anomalies in real-time during the cultivation. In addition, the estimation procedure is robust in the case of noisy observations, and it is not limited to a few phenological stages.

1. Introduction

Knowledge of the current state of an agricultural crop (i.e., its phenological state) is a powerful tool in the context of precision farming. It can be used to predict the future yield production by using agronomical models, to assess irrigation requirements, to plan or trigger fertilisation activities, to detect anomalies due to pests, etc. Remote sensing represents a unique technology to measure vegetation status over large areas with short revisit periods. For this purpose, time series of vegetation indexes (VI) can be used to estimate different dates along the cultivation cycle, like Start of Season (SoS), Peak, and End of Season (EoS) [1,2,3,4].

The Normalised Difference Vegetation Index (NDVI) is one of the most employed VIs in precision farming. With the goal of extracting seasonal stages from the NDVI data, some methods are based on thresholds, which are used to separate different phenological states (greenup, maturity, senescence, etc.). In general, these thresholds are associated with changes in the curvature of the typical NDVI time series. In threshold-based methods, it is common that noise causes incorrect state estimates and false alarms. To minimise the influence of noise (independently of its source), some type of smoothing or fitting can be used [2,5,6,7,8,9,10]. Fitting-based methods employ the parameters of the fitted functions (e.g., maximum or minimum rate of change in curvature) to derive particular phenological stages. However, as exposed in White and Nemani [11], all of these methods cannot work in real-time applications because they require the whole time series to obtain the estimates of dates and stages. Other methods propose a real-time or near real-time operation [11,12,13,14]. For example, in Suwannachatkul et al. [12] the Kalman filter (KF) is proposed as a real-time method to estimate the current state of rice fields, but this method only provides a few particular stages, although the process being monitored is continuous. In White and Nemani [11], the authors present a land surface phenology estimation and conclude that, in the case of crops or homogeneous regions, the use of growth models combined with remote sensing observation provides better results. In our case, we propose exploiting dynamic techniques that combine the information provided by a prediction model and remote sensing data. This approach is employed in many fields, such as robotics [15] or target tracking [16], among others, for optimal estimation of noisy state variables.

In this study, a simple model is used to provide a prediction of the next phenological state. The prediction model takes into account the previous state and the elapsed time between the previous state and the current time. This prediction is then combined with the value inferred by the NDVI data, thus providing a real-time estimation. From a farm manager perspective, it constitutes an interesting product for different uses, especially those requiring critical temporal sampling [17]. The information provided by prediction and observation can be combined by different types of estimators. One of the most used has been the Kalman filter, which obtains the optimal solution [18]. However, it is restricted to linear models and processes described by Gaussian distributions. In order to avoid these limitations, the particle filter (PF) [19] has been selected here to estimate the phenological states.

In order to illustrate the potential of this technique, in this work, we show results, obtained over a set of rice plots, of two different applications: (1) estimation of the phenological state after every acquisition; and (2) estimation of the date of End of Season (EoS).

The text is organised as follows. Section 2 describes the methodology, including the particle filter theory and the particular implementation used here. The available data set and the test site are introduced in Section 3. Results are presented in Section 4 and discussed in Section 5. Finally, Section 6 summarises the main points.

2. Methodology

Systems that evolve over time are considered as dynamic systems and can be characterised by a set of differential equations employed to describe mathematically their physical behaviour. One or more time dependent variables ($x_1(t)$, $x_2(t)$, ..., $x_n(t)$) are direct or indirect characteristics that completely describe the system state at time t. Instead of regarding variables as functions of time, they can be thought as coordinates of points (states) in an n-dimensional domain. The space of all possible points for the system is the state space. One of the most relevant aspects of the dynamical methods is that the future states are inherent to the model description. In other words, the next state of the system at instant $t+1$ could be inferred by the model from the state at instant t. Generally, a dynamical system is described by the following equations:

$$x_{k+1}=f(x_k)+v_k,$$

(1)

$$y_k=h(x_k)+e_k,$$

(2)

where Equation (1) describes the evolution of the state with time, and Equation (2) relates the (noisy) observation to the state. Function $f()$ is known as the dynamic model or prediction model, i.e., a representation of the expected behaviour, linking $x_k$ and $x_{k+1}$. In this notation, $x_k$ is the current state and $x_{k+1}$ is the next state. The observation $y_k$ and the state are related by $h()$, which is known as the observation model. Furthermore, $v_k$ and $e_k$ represent the model error and the measurement error, respectively, which correspond to stochastic noises defined by their known probability density functions (PDF) $p_{v_k}$ and $p_{e_k}$. In this case, functions $f()$ and $h()$ do not depend explicitly on time t, hence resulting in a time-invariant system. This means that we do not have to redefine the expressions of $f()$ and $h()$ for future scenarios.

Therefore, the estimation process is supported by the dynamics and not just by the experimental measurements (observations) themselves. In this context, the phenological development of a particular field can be regarded as a process of a dynamical system, i.e., the crop. The state space will be defined by one or more features that represent the state of the crop. For this work a one-dimensional state space, using phenology as feature, is employed. Therefore, Equation (1) provides the prediction of the next phenological value given the current value, and Equation (2) relates phenology (i.e., state) with the remote sensing observations.

2.1. Particle Filter (PF) Theory

Prediction and observation models are used to infer the system state. Estimation algorithms use the relationship between both sources Equations (1) and (2) to provide a single estimate by their combination. The PF algorithm is selected because it is not constrained to work with linear models and Gaussian distributions.It is based on Monte Carlo methods [19,20] and recursive Bayesian Sequential Estimation [21,22] that provides the posterior PDF ($p(X_{k+1}|Y_{1:k+1})$) at instant $k+1$, combining all the available information (prediction and observations). In order to provide the posterior PDF, the process is separated in two steps: prediction (Equation (3)) and update (Equation (4)).

In the prediction step, the model is used to generate the prior PDF that defines the most likely state at time k without introducing the observation:

$$p(X_{k+1}|Y_{1:k})=\int p(X_{k+1}|X_k)p(X_k|Y_{1:k})dX,$$

(3)

where $p(X_{k+1}|X_k)$ represents the transition probability from state $X_k$ to $X_{k+1}$ at time $k+1$, and $p(X_k|Y_{1:k})$ is the posterior PDF at time k. Finally, when the observation $Y_{k+1}$ is available, the prior probability density function (Equation (3)) is updated via Bayes’ rule:

$$p(X_{k+1}|Y_{1:k+1})=\frac{p(Y_{k+1}|X_{k+1})p(X_{k+1}|Y_{1:k})}{\int p(Y_{k+1}|X_{k+1})p(X_{k+1}|Y_{1:k})d{X}_{k+1}},$$

(4)

where $p(Y_{k+1}|X_{k+1})$ expresses the likelihood function that describes the measurement model, and $\int p(Y_{k+1}|X_{k+1})p(X_{k+1}|Y_{1:k})d{X}_{k+1}$ is a normalisation constant. The state of the system is determined by the most likely value of the posterior PDF. Considering the Markovian assumption, the posterior PDF $p(X_{k+1}|Y_{1:k+1})$ can be obtained recursively from the PDF $p(X_k|Y_{1:k})$ calculated at previous state k. In case the initial distribution PDF $p(X_0|Y_0)=p(x_0)$ is unknown, a uniform distribution over the state space can be considered.

The PF is used to approximate the posterior PDF (Equation (4)) with a set of N samples when it is not analytically tractable. Each sample or particle represents a specific state of the system ($X_{k+1}^i$), and each one has a probability value (named weight, ${\omega }_{k+1}^i$) given by the PDF. In other words, the PDF is represented by a set of N particles and their weights:

$$p(X_{k+1}|Y_{k+1})\approx \sum _{i=1}^N{\omega }_{k+1}^i\delta (X_{k+1}-X_{k+1}^i),$$

(5)

where $p(X_{k+1}|Y_{k+1})$ is the true posterior PDF, $X_{k+1}^i$ is the i-th simulated sample (particle) and ${\omega }_{k+1}^i$ is the weight of i-th simulated sample (particle). N represents the number of particles employed in the simulation, and δ is the Dirac delta function. If N is sufficiently large, Equation (5) approximates the true posterior.

The whole procedure is described in

Table 1. Particle Filter (PF) Sequence.

(1) Initialisation	Generate N samples of $x_0$ from the initial PDF $p(x_0)$.
(2) Prediction	Obtain the sample of $x_k^i$ from the transition PDF $p(x_{k+1}|x_k^i)$.
(3) Measurement step	Compute the likelihood function. $p(y_{k+1}|x_{k+1}^i)$.
(4) Update	Evaluate the importance weights from likelihood function. ${\omega }_{k+1}^i={\omega }_k^{i}p(y_{k+1}|x_{k+1}^i)$.
(5) Normalisation	Normalise the weights ${\omega }^i=\frac{{\omega }_{k+1}^i}{{\sum }_{i=1}^N{\omega }_{k+1}^i}$.
(6) Resampling	The effective number of particles ($N_{eff}$) provides a measure of the number of particles with significant weight representing the posterior PDF. If this number is lower than a provided threshold ($N_{thrd}$) they are redistributed where the PDF is more likely. Reset to ${\omega }^i=1/N$.

and Figure 1. The PDF evolves in time according to the transition model $p(X_{k+1}|X_k)$, which is applied to each particle generating a new distribution named prior PDF. Once the observation is available, the prior PDF is updated to compute the posterior PDF. An important step in the PF is the principle of sequential importance sampling, used for selecting the particle weights [23]. Generally, it is difficult or impossible to draw samples from the posterior PDF $p(X_{k+1}|Y_{1:k+1})$. To avoid this difficulty, importance sampling generates particles $x_{k+1}^i$ from a function $q(x_{k+1}|y_{k+1})$ known as a proposal distribution (or importance density) and assigns the weights (importance weights) according to

$${\omega }_{k+1}^i\propto \frac{p(X_{k+1}^i|Y_{k+1})}{q(X_{k+1}^i|Y_{k+1})},$$

(6)

where $q(X_{k+1}^i|Y_{k+1})$ is the proposal density. As presented in Arulampalam et al. [24], the proposal distribution is factorised such that

$$q(X_{k+1}|Y_{k+1})=q(X_{k+1}|X_k,Y_{k+1})q(X_k|Y_k).$$

(7)

By substituting Equations (4) and (7) into Equation (6), the following equation is obtained

$${\omega }_{k+1}^i\propto {\omega }_k^i\frac{p(Y_{k+1}|X_{k+1}^i)p(X_{k+1}|X_k)}{q(X_{k+1}^i|Y_{k+1})},$$

(8)

and the weight set satisfies:

$$\sum _{i=1}^N{\omega }_{k+1}^i=1.$$

(9)

Equation (8) expresses the basic principle of the sequential importance sampling filter [24]. Finally, the most convenient and frequent is to choose the transition prior as proposal importance density (Equation (10)) to simplify the calculus in Equation (8) so that the weights derive on a recursive set given by Equation (11):

$$q(X_{k+1}|Y_{k+1})=p(X_{k+1}|X_k^i),$$

(10)

$${\omega }_{k+1}^i\propto {\omega }_k^{i}p(Y_{k+1}|X_{k+1}^i).$$

(11)

Table 1 and Figure 1 summarise the key steps to implement the PF algorithm. First, the particles are distributed following the initial PDF $p(x_0)$. Afterwards, the iterative process starts in the prediction step (point 2 in Table 1). The particles that approximate the posterior PDF at instant k go through the prediction model. It gives the new values of the state variables at instant $k+1$ that represent the prior PDF $p(x_{k+1}|x_k)$. In order to compute the likelihood function, the state variables are transformed to the observation domain using the observation model in the measurement step, Equation (2). Given the observation, $y_{k+1}$ and $p_{e_k},$ the likelihood function is computed with the purpose of updating the weights of the particles. Finally, the weights are normalised to satisfy Equation (9). The result is the new posterior PDF from which the estimation is inferred. This value can be determined in different ways, as the most likely value of the PDF (i.e., the mode), the median or the mean, among others. With the iterative process, the particles suffer a degeneration effect: the number of particles with a representative weight tends to be low. A measure of the effective number of particles can be approximated by Equation (12) introduced in Bergman [21]:

$$N_{eff}=\frac{1}{{\sum }_{i=1}^N{w}_i^2},$$

(12)

with a small value indicating a severe degeneracy problem. In order to solve this problem, a technique known as resampling is employed. Different resampling methods have been proposed, as presented in [24,25]. In this work, the residual-systematic resampling (RSR) is implemented [26]. Its goal is duplicating the number of particles in the areas with high normalised weight and discarding the particles with low normalised weight. In this step, N particles are redistributed over the areas where the particles present more weight. Finally, the weights are reset to uniform values ($1/N$). The iterative process continues at point 2.

2.2. Particle Filter Implementation

Provided that agricultural crops are systems that evolve over time, they can be studied in the dynamical framework introduced in the previous subsection. Consequently, filtering techniques can be applied to estimate the phenology of a crop [27,28]. The state space is naturally represented on the Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie (BBCH) scale, which defines the growth states along the life cycle of many cultivated plants [29]. It uses a continuous code range between states 0 and 100, where 0 is associated with sowing and 100 refers to the harvest. In addition, the states are grouped in 10 main intervals or stages, defined by the tens in this scale: germination (0–9), leaf development (10–19), tillering (20–29), stem elongation (30–39), booting (40–49), heading (50–59), flowering (60–69), development of fruit (70–79), ripening (80–89), and senescence (90–99).

This methodology should be applicable to any class of crop in any location defining the corresponding prediction and observation models as shown in Figure 1. In this subsection, an example of the implementation is given for a particular class of rice.

2.2.1. Crop Phenology Model

The first goal consists in finding a model valid for characterising the phenological evolution of rice. Since such a model will be used to provide an estimate of the next state reached by the crop when the current state is known, a mathematical approximation of the expected behaviour of a rice field has to be obtained. As it was mentioned in the Introduction, changes in phenology of rice are mostly driven by temperature evolutions [30], but here we propose a simplified model in which phenological changes are a function of time and not of temperature. Consequently, this model is not the best because temperature information is not taken into account, but it is improved with respect to the ones based only on observations by the combination of prediction and observation during the estimation step. Moreover, this framework would be also appropriate if temperature data were available, since temperature could be included in the prediction model as a control variable.

In order to provide a prediction model that represents the phenological development, a set of ground truth data is employed to fit a function. A detailed description of the ground truth data will be presented in Section 3, but here we use them to illustrate this step. Considering all parcels and all ground acquisitions between 2008 and 2013, we dispose of a total of 656 measures, which are all shown in Figure 2. Due to the fact that the parcels were sowed at different dates, it is impossible to fit the data using the day of year (DoY) as time reference. In order to extract the common behaviour of the crops, the time reference is fixed to the sowing date, and hence the ground data are represented using the day after sowing (DaS) instead of the DoY. Using this representation, all the parcels start at the same time in the phenological state 0 (germination). In the early states, between leaf development and the end of stem elongation phase (states 10–39), phenology increases at a constant rate. Then, at the beginning of the booting phase (state 40) this rate of change increases until plants reach the ripening phase (state 83). Finally, from ripening to senescence phase (state 90), phenology is characterised by a small change rate. Such a temporal evolution can be represented mathematically by a linear trend up to state 40, followed by a sigmoid function, as proposed in the following expression:

$$x(t)=\left\{\begin{array}{cc}mt+n& t<t_c\\ a+\frac{b}{1+exp(-r(t-t_0))}& t\ge t_c\end{array},\right.$$

(13)

where t is time measured in days after sowing; $x(t)$ is the phenological state at time t; m and n are the slope and offset of the line, respectively; $a+b$ is the maximum phenological value, being a the initial background value; $t_0$ is the inflexion point and r the ratio of the logistic curve; $t_c$ is the day at which the function is switched.

In order to fit Equation (13) to the available data, all parameters are tuned simultaneously by minimising the mean square error (MSE). A graphical representation of the fitted function is shown in Figure 2 by a solid line.

Taking derivatives of Equation (13), one can characterise the system by a state-space formulation. After a few algebraic manipulations, the state equation becomes:

$$\dot{x}=\left\{\begin{array}{cc}m& t<t_c,\\ \frac{r}{b}\left(x(t)-a\right)\left(b-x(t)+a\right)& t\ge t_c,\end{array}\right.$$

(14)

To formulate a computational solution, a discrete-time model is required. The discrete version of Equation (14) can be obtained easily by following the Euler method to approximate the derivative, considering small values for $\mathsf{\Delta }t$, yielding:

$$x_{k+1}=\left\{\begin{array}{cc}{x}_k+m\mathsf{\Delta }t& x_k<x_c\\ x_k+r\mathsf{\Delta }t(x_k-a)(b-x_k+a)/b& x_k\ge x_c,\end{array}\right.$$

(15)

where $t_k=k\mathsf{\Delta }t$ in which k takes integer values; $x_k$ is the phenological state at instant $k\mathsf{\Delta }t$; $x_{k+1}$ is the phenological state at the next instant; and $x_c$ is the phenological state at which the function is switched.

In Figure 2, we can also appreciate that there is a distribution of the ground truth data around the model (fitted expression), which is just a representation of the average behaviour. The expected variation around the model is incorporated through $v_k$ in Equation (1). In order to characterise the PDF of the noise model ($p_{v_k}$), the state transition variations are considered. Parcels starting at the same phenological state evolve to the same average state but with small differences. The distributions after each transition are employed to describe the transition noise itself. The expected value of the transition (mean) is obtained applying the evolution model to the phenological state (before the transition occurs). Both values are used to model $p_{v_k}$ with a Gaussian distribution centred in the value given by the model and the variance observed in the transitions.

2.2.2. Observation Model

Besides the availability of the model to make predictions of next states, producing accurate estimates requires the incorporation of observations provided by a sensor, e.g., a satellite acquiring images, in the estimation process. As mentioned in Section 1, NDVI is strongly related with vegetation status, so it can be used as indicative of the phenological state. At this point, the observation model (2) will represent the relationship $h()$ between the NDVI acquisitions, $y_k$, and the state, $x_k$, of the system at a particular instant.

In order to find the expression of $h(x_k)$, all the available acquisitions for a set of rice fields for several years (379 measurements in total) were considered. The NDVI values acquired by the sensor are represented as a function of phenological state in Figure 3. It is observed that, up to the tillering (around state 20), measurements exhibited low NDVI values, followed by a steep increase and saturation around 0.8. Finally, they drop to low values again at states 70–80, approximately.

In some studies, it has been demonstrated that NDVI time series can be approximated by fitting double logistic curves [5,31,32,33] or asymmetric Gaussian model functions [6] to minimise the measurements noise, or noise can be smoothed out using the Savitzky–Golay filter [34]. Here, we will follow the approach based on fitting a double logistic curve, hence:

$$h(x_k)=c+d\left[\frac{1}{1+exp(-r_1(x_k-f_1))}\right.+\left.\frac{1}{1+exp(-r_2(x_k-f_2))}-1\right],$$

(16)

where c is the minimum NDVI value and $c+d$ the maximum NDVI value, $r_1$ and $f_1$ are parameters of the first logistic curve, and $r_2$ and $f_2$ correspond to the second logistic curve. As in the previous case, all parameters are fitted to optimise the MSE.

The observation model (2) is completed by modelling the so-called observation noise $e_k$. We have seen in Figure 3 that a particular phenological state presents a range of NDVI values acquired by the sensor. The observation noise is modelled by a Gaussian distribution with zero mean and variance derived from the data dispersion with respect to the NDVI value given by Equation (15).

2.2.3. Estimation

To estimate phenological states, a uniform prior distribution is taken initially, $p(x_0)$, between states 0 and 50; as in the first observation, the crop will not be in later states. Then, the first NDVI image is incorporated and $p(y_{k+1}|x_{k+1})$ is calculated for every particular case using Equation (16). Next, Equation (11) is used to compute the probability a posteriori. The current state corresponds to the most probable value in the PDF. In the first estimation, the state is determined only with the observation since the prior distribution is uniform. For the next observation, in the prediction step, particles evolve following the transition equation from the instant of the previous observation until the current one. Then, the process is repeated each time an observation is acquired. Regarding the estimation of the date of future states, the prediction model can be used at any time to estimate in which day a specific phenological state will be reached. In that case, additional observations are not incorporated since it is a future event.

In order to employ this methodology over other kind of crops, the next steps have to be followed. First, the state-space must be defined by a set of variables ($x_k$) that characterises the system at any time. In our case, a 1D space is used in which the state variable is phenology, but other physical variables such as leaf area index (LAI), biomass, etc. could be employed to derive a new state-space. The second step involves the definition of the prediction model. In order to provide function $f(x)$, a set of phenological measurement for different parcels and different years are used. Using a fitting strategy, the measured values are adjusted to provide the evolution model and the deviations of the data define the noise process $p_v$. The third step regards the observation sources. The evolution of the system is observed by a particular sensor given a set of measurements. Each state responds with a certain value to the observation. Given a set of measurement for different states of the system, the equation $h(x)$ can be derived. The observation noise is defined considering this perturbation. Once prediction and observation models are defined, the implementation of the particle filter is applied in the same way that is explained in Section 2.1.

3. Data Set and Test Site

This methodology has been tested with data acquired over rice fields located in Sevilla, southwest SW Spain (see Figure 4), during 2008–2013 (2012 is not included). In this area, rice is cultivated from May to October, with just one crop per year. The phenological information acquired during ground campaign measurements have been provided by the local association of rice farmers (Federacion de Arroceros de Sevilla) on a weekly basis over eight to 13 specific parcels, with areas ranging from 3 to 17 ha, depending on the year (see

Table 2. Set of parcels with auxiliary phenological information available.

Year	2008	2009	2010	2011	2013
Number of Parcels	11	13	13	9	8
Images per Parcel	15	16	15	14	6

). In addition, both sowing and harvest dates are known for each parcel. This information was used in the generation of the model and in the validation of the results.

Table 2 shows the total number of NDVI images available per parcel and per year, but only the images during the campaign are used. The images were acquired by Landsat 5 and Landsat 7 with a pixel size of $30\times 30$ m. All the images were already geometrically and radiometrically corrected [35,36]. Some images were discarded due to excessive cloud cover, so the temporal basis (time gap between successive images) is not constant. In the best case, we can obtain an eight-day basis, combining products of both sensors.

To validate the methodology, the set of parcels used to build the model has to be different from the set employed to obtain estimates. As the data set is small, a cross validation method has been employed in this analysis [37]. In this procedure, one case (i.e., all the data from a single parcel) is withheld and models are generated with the remaining cases. Parameters of both Equations (15) and (16) are obtained applying the cross validation method. Therefore, they will be different for each case of study. Then, estimates are computed for the withdrawn parcel. This allows us to maximise the information used in the generation of the model.

4. Results

4.1. Phenological State Estimation

The first case study is focused on the estimation of the current phenological state each time NDVI information is available. For all the examples, the particle filter was implemented with a total of 1000 particles.

The estimation is carried out over the available parcels and images for the five years, resulting in a total of 379 estimates. Results are shown in Figure 5. A high determination factor $(R^2=0.93)$ is obtained, with an root mean square error (RMSE) of $6.6$ states.

In order to illustrate with an example the performance obtained by the proposed dynamic tool, the following evaluation is made for a single parcel (see Figure 6). Phenological states are computed employing two different approaches. On the one hand, the proposed methodology is used to provide the estimation, and the results are shown with the blue line. On the other hand, the estimation is inferred using only observations (without taking into account the model predictions), showing the results with the red line. In the latter case, the states are determined by the maximum values of the likelihood PDF $p(y_k|x_k^i)$ that are derived from Equation (2).

4.2. Prediction of Key Dates

The second case study is focused on the estimation of the date named end of season (EoS), and it also serves as a comparison with other methods widely used in the literature. Regarding the other methods, hereafter called fitting-based methods, the well-known TIMESAT software [38] was used to adjust the available time series in three different ways: with an asymmetric Gaussian curve, with a double logistic function, and by applying the adaptive Savitzky–Golay filtering.

In the fitting-based methods, the estimation of the EoS is computed as the date in which the fitted curve reaches a reference threshold. Instead of directly using the NDVI, these methods were applied to the time series of the $NDV{I}_{ratio}$ defined in Equation (17):

$$NDV{I}_{ratio}=\frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}},$$

(17)

where $NDV{I}_{min}$ and $NDV{I}_{max}$ are the minimum and maximum values of the NDVI time series, respectively. The transformation to $NDV{I}_{ratio}$ allows us to fix a common threshold to estimate the EoS, independent from the specific land cover and geographical location [39]. The EoS was determined as the day in which the $NDV{I}_{ratio}$ reached 50% in a downwards direction.

The correct use of TIMESAT software to provide seasonal estimations employing NDVI time series requires a uniform time sampling of the input data. The presence of clouds forced us to discard some images, so the NDVI was interpolated to provide an eight-day time basis. For some time series approaches, it must be noted that uniform sampling is a limiting factor for solving problems directly in the time domain. In contrast, in a state-space approach, the most relevant information is the behaviour projected in the state space, i.e., the trajectory. For this reason, such an approach is advantageous for scenarios with a non-uniform sampling rate, as it often occurs within optical satellite data. In addition, setting thresholds is not necessary in our approach. The only requirement is knowing the correspondence of the EoS in the phenological scale. According to Meier [29], the EoS corresponds to the phenological state 92.

From the total set of available data, we did not consider in this test year 2013 because only a few images were available and, consequently, fitting-based methods can not provide a useful adjustment to estimate the EoS. Moreover, in the 2011 campaign, only four rice fields reached state 92 during the ground campaign. Therefore, a total of 42 cases (rice parcels monitored during whole campaigns) are used for this test.

In this evaluation, we employ only the first three NDVI images of the campaign to estimate the crop state, and then the prediction model is used directly to provide the EoS estimate. This means that the forecast of EoS is given between 40 and 60 days before it occurs, i.e., it is for more than a month.

Figure 7 shows a comparison of the results obtained by using the four methods: (a) asymmetric Gaussian curves, with a determination factor ($R^2$) of $0.4$; (b) double logistic functions, with $R^2=0.5$; (c) adaptive Savitzky–Golay filtering, with $R^2=0.64$; and (d) dynamical approach based on particle filter, with $R^2=0.77$.

To evaluate the accuracy of all methods for a reduced number of images, the following comparison is done. For each parcel in 2009, the EoS is estimated using TIMESAT software and compared with the results obtained by the dynamical framework. Due to the requirements of a uniform time sampling, the NDVI values are interpolated to provide a total of 38 input data with a constant sampling rate of eight days. The RMSE and the maximum absolute error (MAE) are presented as measurements of the accuracy obtained with the different methods. Results obtained with TIMESAT are summarised in

Table 3. RMSE and MAE obtained for the EoS estimation using TIMESAT.

Images Employed		16	15	14	13	12
asymmetric Gaussian	RMSE	7.2	8.1	9.2	11.0	14.0
MAE	13	21	35	35	43
Double-logistic	RMSE	7.2	8.1	9.0	11.0	14.0
MAE	13	21	37	39	43
Savistzky-Golay	RMSE	10.2	11.0	12.4	14.0	17.7
MAE	20	26	45	47	69

. The first row shows the number of images used in the interpolation step. The results obtained when the whole set of images (16) is taken into account are shown in the first column. The consecutive columns show the results when one, two, three and four images are removed before the interpolation step is applied. The process to estimate the EoS when the number of images decreases is explained next. First, one image is removed randomly from the whole set and the remaining ones are interpolated to provided a constant sample rate, i.e., 38 NDVI input values. The EoS is estimated using the TIMESAT for each method. The process of removing and estimating is repeated a total of 30 times, providing 30 different estimates for each parcel and method. With these values, the RMSE and the MAE are computed. Then, the same approach is employed, removing two, three and four images from the original set.

Table 4. RMSE and MAE obtained for the EoS estimation using a dynamical approach.

Images Employed	3
RMSE	8.3
MAE	24

shows the results obtained using the method proposed in this work. For this case, the interpolation step is not necessary. In order to provide the EoS, two steps are applied: first, three NDVI images are combined with the prediction to produce the estimation of the current phenological state. Second, the prediction model is employed to compute the predicted EoS.

4.3. Estimation over Other Types of Rice

In this last example, we study the reliability of this methodology when a crop grows in a different way from that expected (it grows faster or slower than usual). In such a case, our prediction model can not follow the changes because it does not include any other parameters (external information) except time. To evaluate this effect in the final estimates, we test the method over a crop with a faster phenological evolution (120 days in total) since the model is designed for rice whose cycle is around 150 days. As real data from a faster crop are not available, we simulate it by compressing the time scale.

The proposed methodology is applied in the same way as before, without making any other changes, and results are compared with the estimates of the same crop without changes in the time scale. In Figure 8b, we show the estimation of the phenological states when time between acquisitions is compressed by 20% and in Figure 8a without such a compression. Each time a new observation is available, two estimations are computed: the estimation given by the prediction (green triangle) and the estimation given by the observation (red circle). The particle filter combines both estimations to provide the final phenological state at each step.

The same process is extended to the whole set of parcels providing a total of 379 estimations (see Figure 9). In this case, comparing with the results in Figure 5, the $R^2$ is slightly reduced to $0.90$, whereas the RMSE doubles ($11.2$).

5. Discussion

5.1. State Estimation and Prediction

The estimation of the current phenological state for a specific crop field describes one of the main objectives of this methodology. Hence, the first proposed evaluation is addressed to quantify the ability of the algorithm to estimate such a variable. The importance of providing a real-time estimation has to be remarked here. To be employed during the cultivation campaign, information has to be delivered as soon as possible to the final users so as to conduct the required actions. Although this feature has not been directly evaluated in a final application in precision farming, it represents a unique characteristic of this method and distinguishes it from alternative state-of-the-art techniques.

To test the quantitative performance, the state estimations have been compared with the reference data provided by the on site campaign (see Figure 5). In general, the estimation accuracy improves when more observations are available. However, the performance is worse from state 40 to 60 than for the rest. This is due to the lack of sensitivity of the observations (i.e., NDVI) to phenological changes in this interval (see Figure 3), so the estimation is mainly determined by prediction. It should be emphasised that observations alone would not be able to provide any different estimate in this interval, but the proposed methodology can cope with these conditions. Thanks to the model and the previous estimation, we can keep a good tracking between these states (see results between DoY 220 and 260). Finally, after state 60, the NDVI provides useful information again for updating the prediction and improving the performance of the estimation. The results presented in Figure 5 demonstrate that this methodology provides reliable information about the crop status all along the cultivation campaign.

Moreover, a comparison between the PF estimation and the direct inversion of the observation proves the ability of the approach to solve NDVI evolution ambiguities. An interesting feature of these results is the impossibility of estimating some phenological states using only NDVI observations. As it was mentioned, NDVI exhibits ambiguity between both initial and final states (both of them are characterised by low values), whereas the middle states present high constant NDVI values. Consequently, the third estimation in Figure 6 is totally incorrect. The likelihood PDF of an NDVI value around 0.5 (which is present both at the initial increase and the final decrease) is illustrated in Figure 10. The resulting bimodal distribution justifies the erroneous estimates provided by the observations alone. On the contrary, the particle filter can deal with bimodal distributions exploiting the prediction step to weigh with low values all improbable transitions, hence avoiding these unrealistic jumps. This is the main reason of using the particle filter instead of an extended Kalman filter (EKF) which only works with Gaussian distributions.

An alternative, and useful in practice, output product from the algorithm is the estimated date of arrival or prediction of key dates. These products are based on a reduced set of observations and the evolution model to infer at which date a specific state will be reached. A comparison in the precision provided by fitting-based methods and the proposed methodology has been conducted. The TIMESAT software has been selected to this aim providing predictions for the date of arrival using several methods, all of them based on curve fitting over the NDVI time series. The comparison strategy consisted of estimating the date of the same state, in this case the end of season (EoS) was selected, with all methods.

This comparison shows the evident contribution of the proposed methodology, which outperforms the other methods and is able to provide a good estimate of the EoS (RMSE = 8 days) with more than 40 days of anticipation. The accuracy of the fitting-based methods is strongly related with the number of images employed, providing better results with a large number of data. The adaptive Savitzky–Golay filtering presents the best results when all years are taken into account. The exception for this case is 2009. This is because some parcels exhibited an increment of the NDVI during some months before the sowing date. Probably, it is an effect of natural vegetation growing in the location of the parcel that is removed before the sowing. The asymmetric Gaussian and double-logistic fitting are less sensitive to these perturbations. The method proposed here makes use of the prediction model to reduce the number of observations employed and the interpolation step is not necessary. In this case, using just three images, results are similar to the case of employing 15 images with the fitting-based methods (see Table 4).

5.2. Methodology Generalisation

This subsection provides a brief analysis on the generalisation process and the ways for it to be applied. The presented applications were tested using the same set of crop fields both at the modelling and evaluation steps. However, the model is designed to be employed over alternative or future campaigns. Due to internal (type of rice) or external (temperature, soil characteristics, pests, etc.) conditions, the phenological development can be faster or slower with respect to their expected behaviour, i.e., the evolution model. For this reason, it is mandatory to evaluate how these factors affect the estimation process.

To perform such an analysis, a simulation of a faster development has been derived from the original set. Figure 8 shows the comparison between the original phenological evolution (150-days evolution) and the faster development (120-days evolution). In this case, there exists a general trend, caused by the evolution model, providing underestimated states. The prediction gives lower phenological values due to the slower evolution of the model in contrast to the simulated 120-day rice. This delay is partially compensated in the prediction with the filter when incorporating observations. As expected, the results show that when using the wrong model, the estimations are worse. The deterioration is due to the low sensitivity of the NDVI between states 40 and 80 combined with a bad prediction. This specific range of states accumulates most of the error based on the prediction since the observation does not provide enough information to update the state correctly. However, we could use another model for crops that evolve faster, and hence improve the prediction step. Another alternative is to incorporate a control variable in the model that contains information about the velocity of the development of the crop to lead the prediction, for instance the accumulated temperature.

5.3. Summary of Advantages

The dynamical approach proposed in this work enables the use of a simple prediction model to estimate the phenological state. This model takes into account the previous state and the elapsed time between the previous state and the current time. It works properly when it is applied under similar conditions to those under which it was created (same crop, same location, etc.). This methodology has substantial advantages over conventional methods based on fitting or interpolation. The first one is real-time evaluation, i.e., the estimation is provided each time an NDVI image is received, instead of waiting until the whole set of images is available. Moreover, the estimation is not just based on observation data, but it is optimally combined with the estimation predicted by the model. This derives into the second advantage: supplying a more robust system under persistently noisy observations if it is properly defined. The dynamics of the process can keep a satisfactory estimation performance even when the observation is highly noisy. The third advantage regards the continuity of the output. The model defines a wide set of phenological states instead of just a few of them. The fourth advantage is that the methodology can be used to provide future estimations. This opens a new set of different applications allowing the prediction of the day at which a crop reaches a specific state.

5.4. Perspectives and Future Research Lines

The potential of a dynamical approach based on the PF to provide real-time phenological states and forecasts of state arrival dates has been introduced in this manuscript. In a future work, it would be interesting to test the methodology over a bigger set of data and evaluate the scope of its generalisation to other types of crops. Moreover, the proposed technique can be used in a more global scale, for example in ecological studies. In this domain, there has been a significant growth of interest in phenology in the last few years due to shifts in the timing of different phenological phases in plants connected to climate change [40,41,42,43]. Due to the fact that phenology is affected by climatic conditions, it would be interesting to derive information about climate changes using the state estimated by the proposed methodology in a longer temporal study.

Currently, ongoing work is addressed to include complementary sources of information to combine with NDVI, as images acquired with radar satellites, in order to improve the estimations and reduce the time between observations. Data fusion would benefit also from the diversity in sensitivity to the scene features. Moreover, air temperature has been found to be a dominant factor controlling the timing of flowering and other phenological phases [44,45], so the incorporation of control variables in the model (as the accumulated temperature, expressed in Cumulative Growing Degree Day, CGDD) is being studied to improve the prediction step.

The flexibility of the methodology allows us to consider alternatives for the prediction and observation model. For instance, it can be interesting to evaluate the incorporation of a biological growth model to provide the prediction stage under the same estimation context. Although the evaluation of these sorts of models involves the consideration of a large number of parameters and variables, the reliance on the output product could be significantly improved.

6. Conclusions

Dynamic tools have been demonstrated to be useful for providing accurate estimates of the phenological state of rice crops when time series of NDVI data are available. Thanks to the incorporation of a prediction model and an adequate combination of predictions and observations, carried out with a particle filter, this method outperforms other known methods in several aspects. First, the range of possible states to be estimated is made continuous, and not just a few particular states of the cultivation cycles. Second, estimates are provided in real-time, as observations arrive, or even with anticipation. This is in contrast to common methods, based on curve fitting and thresholds, which need to be employed at the end of the cycle.

In order to illustrate the performance of the proposed methodology, estimates along the whole phenological cycle have been obtained for a set of rice fields with NDVI satellite data acquired for five years. A high correlation ($R^2=0.93$, with $n=379$) and low RMSE ($6.6$ states) is achieved. Moreover, estimates of the EoS have been produced with much better results ($R^2=0.77$, $n=42$, RMSE = 8 days) than common methods (provided by the TIMESAT software) and using just three images acquired at the beginning of the campaigns (i.e., with a 40–60 days anticipation). In addition, we have shown that despite the simplicity of the prediction model, thanks to the methodology, we are able to provide good estimations when anomalous behaviours (i.e., not covered by the prediction model) appear.