The dynamics of wetland cover change using a state estimation technique applied to time-series remote sensing imagery

ABSTRACT Monitoring the dynamics of inundation areas in wetlands over contiguous years is important because it influences wetland ecosystem monitoring. However, because the variable nature of wetlands tends to hamper monitoring change analyses, the potential for misinterpretation increases. The Kalman filter (KF) or extended Kalman filter (EKF), which uses recursive processing based on the former information, can be applied to time-series remote sensing imagery. In the experiment, a periodic triangle function of two modulated parameters is treated as the system model, and Normalized Difference Vegetation Index (NDVI) time-series data are used for the measurement model in the correction processes of the state estimation. A decision metric is computed from the mean and amplitude sequence, which results from the state estimation filter. Consequently, an optimal threshold is calculated using a minimum error thresholding algorithm based on a pre-labelled sample. NDVI time-series data from Poyang Lake, China – derived from 250-m Moderate Resolution Imaging Spectroradiometer satellite data obtained from January 2009 to December 2013 – are applied to monitor the dynamics of inundation changes. The results show that the EKF achieves satisfactory results, with 85.52% accuracy in the year 2009, while the KF has an accuracy of 84.16% during that same time.


Introduction
The wetland ecosystem of Poyang Lake, China, has experienced dramatic changes over the past several decades as a result of climate variability and human activity (Feng et al. 2012;Han et al. 2015), especially the operation of the Three Gorges Dam and sand mining (Gao et al. 2014). These changes have had a substantial impact not only on the wetland ecosystem but on other systems as well, including landscape changes, increased flood frequency (Shankman and Liang 2003) and the transmission of infectious diseases to other areas by migratory birds that are sensitive to local climate change (Wu, Lv, et al. 2014). Several studies have been performed to obtain a greater understanding of the variation in wetland cover, including studies that monitor the dynamics of water-level changes via the water extraction method (Hu et al. 2015) and Ice, Cloud, and Elevation/Geoscience Laser Altimeter System (ICESat/GLAS) data (Wang et al. 2013), including CONTACT Chunxiang Cao caocx@radi.ac.cn observations of the bottom topography of wetland areas (Feng et al. 2011). A water-quality measurement retrieval model for wetland areas based on suspended particulate matter concentration parameters and satellite images has also been developed (Wu et al. 2013;Wu, Liu, et al. 2014).
Although dynamic change monitoring based on a single image provides fruitful results in certain applications (Yuan et al. 2015), the ecosystem of Poyang Lake is dynamic; thus, estimates of inundation areas based on individual images acquired on a single date are not representative of the longterm changes that occur in these areas Sun et al. 2014). Moderate Resolution Imaging Spectroradiometer (MODIS) time-series images offer significant advantages because of their wide coverage and ability to monitor wetland areas . Changes in wetland vegetation can be observed using many indexes, including the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (Lunetta et al. 2006;Lhermitte et al. 2008;Hu et al. 2015), wherein long-term NDVI information is obtained from MODIS time-series data. This approach contributes to the objective of monitoring the dynamic changes in inundation areas. Because the fluctuations in NDVI time-series data in wetland areas are due to many factors (e.g. turbidity and water-level oscillations), obtaining the desired results from monitoring analysis is challenging.
Recently, the state estimate technique has been applied in the study of remote sensing (Samain et al. 2008;Kleynhans et al. 2010;Kleynhans et al. 2011;Insom et al. 2015). The Kalman filter (KF), a state estimation technique, has been applied to linear models used in remote sensing for estimating the parameters of dynamic vegetation models for land-cover classes (Samain et al. 2008). In addition, the extended Kalman filter (EKF), which is the nonlinear version of the KF, linearizes the last state estimate of the mean and covariance and is used for land-cover classification (Kleynhans et al. 2010) and change detection of settlement areas (Kleynhans et al. 2011).
In the study area of Poyang Lake, both the KF and EKF have become interesting alternatives for addressing the wetland monitoring system because system models of wetland ecosystems are not well known, which can lead to inadequate results and inaccurate estimates. In this study, we present a technical method that uses NDVI time-series data and state estimation techniques (KF and EKF) to monitor long-term changes in wetland areas over a five-year period (January 2009-December 2013).

Study area
Poyang Lake (28 22 0 -29 45 0 N and 115 47 0 -116 45 0 E), the largest freshwater lake in China (Figure 1), is located in northern Jiangxi Province. The geography of Poyang Lake is divided by Songmen Mountain. The northern portion is narrow and deep and is connected to the Yangtze River, and the southern portion, which is the primary body of water, is large and shallow (Dronova et al. 2011;Feng et al. 2011;. The Poyang Lake region experiences seasonal variation. From April to September, the water area is relatively large, and all of the floodplains are inundated with water from five rivers (Xiushui, Ganjing, Fuhe, Xinjing, and Raohe). The water-level peaks between July and September, at which time the water area reaches >3000 km 2 , supported by backflow from the Yangtze River. Between October and March, the water-level drops, and the inundation area of the lake can decrease to <1000 km 2 (Feng et al. 2012;. Frequent natural disasters (flooding and drought) have had a severe impact on the local ecosystem of Poyang Lake over the past few decades. For example, a serious drought in 2011  significantly decreased the volume of physical water bodies in the study area, resulting in a number of ecological problems, including changes in animal and human habitats.
Another external factor that affects the Poyang Lake ecological system is human activity, such as the construction of the Three Gorges Dam and sand mining. These activities have also caused substantial changes to the wetland ecosystem, which have had an effect on migratory birds (for example ). Therefore, the study of the dynamic changes in wetlands is increasingly important for improving our understanding of wetland ecosystem changes and for enhancing public awareness.

Satellite data
The NDVI time-series data were derived from MODIS 16-day (MOD13Q1 product) images with a 250-m spatial resolution for Poyang Lake and values ranging between 2000 and 10,000 and with a scale factor of 0.0001. The NDVI time-series data were acquired for the January 2009-December 2013 period, constituting a total of 115 time points.
The nearest neighbour algorithm was used to resample five Landsat ETM+ images from the fiveyear period (21 December 2009, 8 December 2010, 11 December 2011 November 2012, and 1 January 2014), which have a 30-m spatial resolution, to match the resolution of the MODIS data. The nearest neighbour algorithm was then applied to assign the pre-defined class label for the calculation of the minimum error thresholding algorithm. Only one Landsat ETM+ image was employed in each study year because the aim of this study is to summarize the proportion of each type of wetland land cover in each year. Thus, one Landsat ETM+ image from the end of the focused year was selected. In addition, to validate the results from the MODIS image analysis of the focus year (2009), two resampled Landsat TM images (12 February 2009 and 26 October 2009) from different seasons were used to select the validation points. These two images were applied because a selection of validating points should absolutely define the class. Therefore, the validation points were selected from these two images, and the selected points did not change during the year.

State estimation technique
Two types of state estimation, the KF and the EKF, are implemented to compare the performance of linear and nonlinear versions with respect to the wetland ecosystem. The proposed method has three major steps, which are described in the flow chart of the research procedure in Figure 2. KF and EKF are applied to obtain the mean and amplitude sequence of the NDVI time-series data. The decision metrics are then computed. A given pixel is classified as a particular class by the minimum error thresholding algorithm based on the pre-labelled pixel.

Kalman filter
The KF (Ribeiro 2004) is typically implemented to estimate the state of a dynamic system. The system and measurement models of a simple linear discrete-time KF are illustrated in the following equations: where x k 2 R n is the system state vector, w k 2 R n is the vector that conveys the system error sources, z k 2 R r is the observation vector, and v k 2 R r is the vector that represents the measurement error sources. A k is the state transition matrix, which reflects the mathematical and physical relationship between system state x k and x kÀ1 . H k is the measurement function, which depicts the relationship between measurement z k and system state x k . The vectors w k and v k are independent sequences of white noise with zero means. Gaussian noise with a zero mean can be described as follows: Additionally, the joint covariance matrix can be described as follows: where the non-negative definite matrix Q is the system noise covariance matrix, the positive definite matrix R is the measurement noise covariance matrix, E ½ is expectation, and the subscript 'T' denotes the matrix transpose. The aim of the KF is to estimate the actual value of x k in Equation (1). Based on the model equation, the five important key steps of the discrete-time KF are summarized as follows: Filtering: wherexðk j kÞ in the above equation is the estimation value of the system state x k , Pðk j kÞ is the error covariance matrix, which quantifies the uncertainty of the estimate defined by E½ðxðkÞ Àxðk j kÞÞðxðkÞ Àxðk j kÞÞ T , and the weighting matrix KðkÞ is the Kalman gain matrix. The KF begins by working with an initial condition valuex 0 and P 0 . The prediction steps of the KF are represented in Equations (5) and (6), which generate an a-priori estimation of the system state at time k. Equations (7)-(9) are the filtering processes or the measurement update equations of the algorithm. They incorporate the measurement value z k into an a-priori estimation to obtain an improved a posteriori estimation, which is the output of the KF at time k.

Extended Kalman filter
To address the filtering problem in case the system dynamics (state and observations) are nonlinear, the system and measurement models can be represented as Because of the nonlinear dynamics, these probability density functions (PDFs) are non-Gaussian. To evaluate the first and second moments, the optimal nonlinear filter has to propagate the entire PDF, which has a heavy computation burden. Rather than propagating the non-Gaussian PDF, the EKF (Ribeiro 2004) considers, at each cycle, a linearization of the nonlinear dynamics of Equations (10) and (11) around the last consecutive predicted and filtered estimates of the state, and for the linearized dynamics, it applies the KF. One cycle of the EKF consists of the following steps: (1) Consider the last filtered state estimatexðk j kÞ.
(5) Apply the filtering cycle of the KF to the linearized observation dynamics, yieldingxðk þ 1 j k þ 1Þ and Pðk þ 1 j k þ 1Þ.
Let FðkÞ ¼ rf k jx ðk j kÞ be the Jacobian matrices of f ð ¢Þ and Hðk þ 1Þ ¼ rh jx ðkþ1 j kÞ denote the Jacobian matrices of hð ¢Þ. The EKF can be stated as follows: Filtering: x It is important to note that the matrices FðkÞ and HðkÞ depend on the previous state estimates; therefore, the filter gain KðkÞ and the matrices Pðk j kÞ and Pðk þ 1 j kÞ cannot be computed off-line as occurs in the KF. Moreover, the matrices Pðk j kÞ and Pðk þ 1 j kÞ do not represent the true covariance of the state estimates.

Initialized state estimation
The KF and the EKF are used to obtain the mean and amplitude sequences. The general filtering problem is formulated using the following model equations: where w k and v k are vectors that represent the system and measurement error noise, respectively, and are sequences of white Gaussian noise with a zero mean. The system model is defined by f k , and h k is the measurement function applied to the prior state vector to obtain the current state estimation.
2.3.3.1. System model. The NDVI sequence for a given pixel was modulated by two parameters. Let x ij k ¼ ½m k a k T denote a system state vector at time k consisting of mean (m) and amplitude (a), and Equation (17) particularizes to Equation (19): where ij denotes the row and column of the image, w k is the system error vector, and A k is the state transition model, which is used for the previous state k to generate the state vector at time k þ 1. The state transition model in this experiment is generated using the periodic signal of a triangle wave represented in terms of a sinusoid equation with period p and amplitude a.
2.3.3.2. Measurement model. The NDVI value of a given pixel at time k denotes the measurement value of the state k that is employed to update the state vector. Therefore, the measurement model in Equation (18) can be written as 2.3.3.3. Initialized parameters. Initialized parameters and observation noise tuning can be applied in many techniques, such as computing the mean from the NDVI sequence of a given pixel for the initial mean value (Chu et al. 2008) and using a trial-and-error method to adjust the standard deviation of the observation model noise. In this study, the initial state parameters (m 0 ; a 0 ) and standard deviation of the observation noise (s v ) were calculated using the fast Fourier transform (FFT) mean and annual components of the sampling data in the study area, which were set to 5% of the NDVI time-series data (Kleynhans et al. 2010). Figure 3 shows the signal energy when the initial state parameters were set by the FFT components 1 and 5.

Decision metric calculation
The decision metric (d), which is determined by the summation of a mean decision metric (d m ) and an amplitude decision metric (d a ), can be computed using the following equation: The mean decision metric (d m ) and the amplitude decision metric (d a ) are shown as the following equations: where m ij and a ij are an average of the mean sequence and an average of the amplitude sequence of the studied pixels, expressed as Equations (25) and (26). In addition, the average of the mean sequence and the average of the amplitude sequence of the neighbour pixel are denoted by m ij d and a ij d , respectively, and listed in Table 1: where m ij k is the mean value of a studied pixel at time k. Similarly, a ij k is the amplitude value of a studied pixel at time k.

Minimum error thresholding algorithm
Generally, the selection of a suitable threshold is a difficult task (Moser et al. 2007). Thus, this study defines the threshold in the probability domain to avoid the influence of robust noise and choice of  interpretation (Zhao et al. 2011). A minimum error thresholding algorithm has been applied in many studies (Yifang and Yousif 2012), including the present one, in which the PDF of a normal distribution was implemented, as shown in the following equation: f ðx j m; sÞ ¼ 1 s ffiffiffiffiffi ffi 2p p e À ðxÀmÞ 2 2s 2 (27) where x represents the sample pixels for each class. In each year, sample pixels are labelled as four types of classes (water, mudflat, submerged vegetation, and emergent vegetation) in accordance with the Landsat ETM+ images, which have better spatial resolution than the MODIS time-series images, which are implemented by a nearest re-sampling algorithm and then layer stacked into the decision metric (d). The coordinates of the sample pixels are mapped to the decision metric (d) to obtain the sampled values. Consequently, the minimum error thresholding algorithm is employed to label the class that each pixel in the decision metric (d) belongs to, as presented in the research procedure in Figure 2.

State estimation results
Using the same system model and the same measurement model, both the KF and the EKF have been evaluated based on root mean square error (RMSE): where y i is the predicted value, which is equal to the KF or EKF estimation value. The termŷ i is the observation value, which is equal to the original NDVI in this experiment, and n is number of data.
If a periodic signal of a triangle wave is used as a system model and the original NDVI value is applied as the measurement value for one studied pixel, the KF and the EKF differ in predicted performance. Four pixels were selected as examples for validating the performance of the KF method and the EKF method. Two of these four pixels are no-change pixels in the water and emergent vegetation classes. The other two pixels are changed pixels from the water class to the emergent vegetation class and from the emergent vegetation class to the water class. Four pixels that correspond to these conditions were selected. This study performed a validation based on the water class and the emergent vegetation class because they are distinct in terms of NDVI values and more clearly illustrate the performances of the KF and EKF methods.
The RMSE of the EKF method is better than that of the KF method in every test case, as presented in Figure 4. Figure 4(a,b) shows the RMSE values of the KF and EKF methods, respectively, for the no-change pixel in the water category. Figure 4(c,d) shows the predictive performance of the KF and EKF methods, respectively, for the pixel that changed from the water class to the emergent vegetation class. Figure 4(e,f) presents the RMSE values for the pixel that changed from the emergent vegetation class to the water class. Finally, Figure 4(g,h) presents the RMSE values of the KF and EKF methods, respectively, for the no-change pixel in the emergent vegetation category. Although they exhibit only slight differences in RMSE values, this difference is important to the overall class separation performance. The EKF should produce better performance than a linear filter or the KF, in particular when monitoring the land cover of a high-fluctuation area, which is characterized by a high incidence of heavy-tail noise.

Minimum error thresholding algorithm results
Supported by sample pixels from Landsat ETM+ images during the end of the five study years (2009)(2010)(2011)(2012)(2013), the number of samples for each class in every year is presented in Table 2. The differences in the number of samples for each class have no effect on the calculation of the minimum error thresholding because the PDF of a normal distribution does not depend on the number of samples, as explained in Equation (27). The coordinate of every sample pixel is mapped into the decision metric (d), and the value from those positions is then used to calculate the threshold of the decision metric. Figure 5 demonstrates the optimizing threshold of the decision metric based on the EKF using a PDF of the normal distribution with data from 2010. The optimal threshold values can clearly where th 1 , th 2 , and th 3 are the thresholds between water and mudflat, mudflat and submerged vegetation, and submerged vegetation and emergent vegetation, respectively; the values are shown in Table 3.

Accuracy assessment
A total of 500 points are randomly generated in the study area with ArcGIS 9.3 software; then, they are imported to a layer-stacked image of the two Landsat TM images from different seasons in 2009(12 February 2009and 26 October 2009. By combining the TM image in RGB and then using 442 visually interpreted, selected, and validated samples ( Table 4) from points that were homogenous within 3 Â 3 pixels in the TM image around those samples, these validated samples are used to calculate the confusion metric. As mentioned before, two Landsat TM images were applied because a selection of validating points should absolutely define the class. Thus, the validation points were selected from these two images to avoid changes during the year. The confusion metrics using the KF and the EKF for the wetland cover classification are presented in Tables 5 and 6, respectively.
The results indicate that the overall accuracy of both methods is greater than 84%: the class separation by the KF produced an accuracy of 84.16% (Kappa coefficient = 0.7895), while the EKF had an overall accuracy of 85.52% (Kappa coefficient = 0.8073). The overall accuracy of the EKF method is higher than that of the KF method; therefore, the fluctuation and oscillation in wetland areas can be overcome by nonlinear techniques. The overall accuracy of the KF and EKF methods agrees with   the state estimation results, but the prediction performance of the EKF method is preferable to that of the KF method in every case. In addition, the user classification accuracy of the KF and the EKF techniques is at least 80.00%, which is usually appreciated for implementing MODIS time-series imagery. The EKF applied to wetland cover classification offers better user accuracy than that of the KF in the classes of submerged vegetation and emergent vegetation. Similarly, the producer accuracy of the EKF in the class of submerged vegetation is higher than that of the KF, thus emphasizing that the EKF, or nonlinear filter, works well with this system. Although the results of the two filtering methods do not differ substantially, this type of precision is beneficial to certain applications, such as wetland management and mitigation planning for natural disasters.
The land-cover changes in Poyang Lake have been substantial and have impacted the wetland ecosystem. The classification results of the annual NDVI time-series data from 2009 to 2013 based on the application of the KF and the EKF are shown in Figures 6 and 7, respectively. The figures depict the proportion of each type of wetland cover as a percentage for each study year. Since the long-term lake water level reached the lowest water levels during 2001-2010 (Zhang et al. 2014), the proportion of the water class for both the KF and the EKF methods has been a small percentage. Then, in 2011, Poyang Lake experienced a serious drought ); according to the classification results, the proportion of water dropped that year. Drought continued to affect the area through the year 2011 and into early 2012. This low lake level was also related to the practice of storing water behind the Three Gorges Dam (Wikipedia 2015). However, Poyang Lake experienced heavy rains during May 2012, which increased the water level in the lake and the corresponding classification results for the year 2012; thus, the proportions of water and mudflat were larger than in the previous year. In 2013, Poyang Lake again experienced a low-water situation (Davison 2014). The classification results from both the KF and EKF methods agree that, following this event, the proportion of water proportion decreased and the proportion of emergent vegetation increased.

Conclusions
This experiment offers an alternative approach to classifying wetland landscapes. It employs MODIS time-series imagery derived from the MOD13Q1 product, which is provided every 16 days, in the study area of Poyang Lake, China. The state estimation technique (KF and EKF methods), which has been implemented in many applications, is modified for the classification of remote sensing imagery. The overall accuracies of the linear (KF) and nonlinear (EKF) filtering differ slightly, and these differences may be an issue of concern with respect to wetland ecosystems in certain applications. Therefore, further research is needed to improve the accuracy of wetland cover classification using state estimation.
Landscape changes are considered to be important for the wetland ecosystem in the study area. The classification performances of the proposed techniques are favourable for this system and are in accordance with past events. For instance, a serious drought occurred in 2011, resulting in a dramatic shrinkage of the bodies of water. Heavy rain in May 2012 caused the proportion of water to rise. A drought in 2013 decreased the area of water.