Application of Multi-periodic Harmonic Model for Classification of Multi-temporal Satellite Data: MODIS and GOCI Imagery

Abstract

A multi-temporal approach using remotely sensed time series data obtained over multiple years is a very useful method for monitoring land covers and land-cover changes. While spectral-based methods at any particular time limits the application utility due to instability of the quality of data obtained at that time, the approach based on the temporal profile can produce more accurate results since data is analyzed from a long-term perspective rather than on one point in time. In this study, a multi-temporal approach applying a multi-periodic harmonic model is proposed for classification of remotely sensed data. A harmonic model characterizes the seasonal variation of a time series by four parameters: average level, frequency, phase, and amplitude. The availability of high-quality data is very important for multi-temporal analysis.An satellite image usually have many unobserved data and bad-quality data due to the influence of observation environment and sensing system, which impede the analysis and might possibly produce inaccurate results. Harmonic analysis is also very useful for real-time data reconstruction. Multi-periodic harmonic model is applied to the reconstructed data to classify land covers and monitor land-cover change by tracking the temporal profiles. The proposed method is tested with the MODIS and GOCI NDVI time series over the Korean Peninsula for 5 years from 2012 to 2016. The results show that the multi-periodic harmonic model has a great potential for classification of land-cover types and monitoring of land-cover changes through characterizing annual temporal dynamics.

1. Introduction

Land-cover information is widely utilized for a variety of purposes in urban and regional planning, natural resources conservation and management, environmental applications, etc. Remote sensing plays a vital role in providing quantitative and qualitative land-cover information as a cost-effective means through image processing techniques.Various methods to extract and interpret valuable information from massive satellite images have been developed for several decades. Land-cover classification is one of the powerful techniques to extract earth surface information and the resulting class-cover types are the basisfor many socio-economic and environmental applications.

Classification is carried out through identifying differentspectral reflectance and remittance properties of different classtypes on the earth’ssurface, resulting in a map of meaningful land-cover themes. A large number of satellite image classification methods and techniques have been introduced and utilized, each having its own advantages and limitations.The research works on reviewing and comparing satellite image classification methods and techniques are performed to guide or select a suitable classification method for a specific study (Lu and Wang, 2007; Li et al., 2014; Jensen, 2009). An appropriate classification method is to be selected based on the requirements of the application. Spectral-based classifiers have been conventionally used in classification of remotelysensed imagery (Zenzo et al., 1987; Mohn et al., 1987). Signal variability over time often results in artifacts in land-covertypesthat cause classification errors. Due to the difficulty in discriminating surface types based on the spectralsignatures at a specific time,multi-temporal techniques characterizing temporal profiles have been proposed an alternative approach. However, relatively, a few multi-temporal features have been explored. Statistical approaches such as principal component analysis(PCA) and seasonal autoregressive integrated moving average (Lovell and Graetz, 2001; Sakamoto et al., 2005) have been designed for the analysis of remotely-sensed time series data. The physical process of multi-temporal data generally shows the characteristics of a time-series phenomenon with seasonal trends (Brims, 1991; Bloomfield, 1976; Jakubauskas et al., 2001). Forthisreason, the temporal patterns ofland coverthrough temporal profile is better analysed in the time domain than in the spatial and spectral domain (Eastman and Fulk, 1993; Samson, 1993; Reed et al., 2006).

Atransformation of multi-temporal data through the determination of the area-under-the-curve or the integral of the temporal profile curve is one of the approaches in time domain (Tucker et al., 1990). Temporal trends in the processes are summarized by a feature that is associated with the temporal accumulation measure. In this study, multi-temporal analysis based on a harmonic model is proposed. The purpose of this study is to provide an efficient methodology forinformation extraction froma sequence of remotely sensed data including data reconstruction and land-cover classification.

Harmonic analysis has been employed in many applications such as seasonal changes for natural and agricultural land use/cover, the phenological patterns of diverse biomes over a long-term period, and characterization of land use/cover (Doran, 1986; Bradley et al., 2007; Jakubauskas et al., 2001). In harmonic analysis, a complex time series can be transformed into sum of a series of sinusoidal waves (terms) if it represents a periodic phenomenon (Davis, 1986).The harmonic components ofthe time seriesfor an individual pixel represent the characterization of temporal variability of local region corresponding to the pixel. Harmonic model characterizes the seasonal variation of a time series by four parameters: level, frequency, amplitude and phase. The level represents the average ofthe overall greenness overthe period and periodicity of response is described by the frequency. The amplitude value reflects the range (size of the wave) and the phase represents the initiation time of variation (offset of the wave) (Lee, 2008).

Multi-temporal remotely-sensed data have been shown to be successful in monitoring seasonal trends in phonological processes. Analyses of the relations between spectral measurements and vegetation related phenomena have been successful as a global vegetation observatory (Lunetta et al., 2006; Pettorellia et al., 2005; Reed, 2006). Especially, reflectance data from Moderate Resolution Imaging Spectroradiometer (MODIS) that is deployed on the Terra and Aqua satellites are obtainable globally on a daily basis and have been shown to have considerable potential for land vegetation studies at a regional and global spatial scale. The GOCI, launched in 2010, has a 2500 km × 2500 km field of view and eight multispectral bands, including visible and near-infrared bands. The GOCI hasspatialresolution of 500 m and temporalresolution of eight-time a day (Ryu et al., 2012). The GOCI also has a considerable potential to understand land covers at reasonable spatial resolution and high temporal resolution and to be utilized forland vegetation studies.

Multi-spectral reflectance data have been transformed and combined into various vegetation indices to minimize the variability due to external factors(Tarpley et al., 1984).The normalized difference vegetation index (NDVI) is the most commonly used vegetation index (VI) and the NDVI versus time profile represents seasonal development history of vegetation. In this paper, the temporal profile of MODIS and GOCI NDVI processesisinvestigated and characterized as a means of reducing dimensionality for multi-temporal classification. The seasonality of vegetation types can be represented with a multiperiodic harmonic model. The parameterization provides physically interpretable valuesto characterize the seasonal development with. Here, using the estimatesthat are obtained through spectral analysis of a sequence of imagery, seasonal periodicity can be incorporated into multi-temporal classification. The resulting classification based on these components reflects different sources of temporal variation.

The availability of high-quality data is also an important issue in the analysis of the time series data (Chappell et al., 2001; Lu et al., 2007; Roerink et al., 2000). In the case of harmonic analysis, noisy data affect the parameter estimation of harmonic components and might possibly produce inaccurate results. MODIS and GOCI data inevitably contain disturbances caused by cloud presence, atmospheric variability and aerosol scattering, and instrument problems (Gu et al., 2009; Huete et al., 2002;Moody et al., 2005;Xiao et al., 2003), which impede further analysis. The pre-processing is required in orderto reconstruct a complete MODIS and GOCI time series. Many pre-processing techniques to reconstruct high-quality data streamhave been proposed for radiometric and atmospheric normalization and correction (Jonsson and Eklundh, 2002; Sakamoto et al., 2005; Ma and Veroustraete, 2006; Reed, 2006). In this study, dynamic compositing method is employed to reconstruct a sequence ofimages having the physical processes with seasonal periodicity. For this, the harmonic model to track seasonal variation through time is used and a Gibbsrandom field (GRF)is used to representthe spatial dependency (Lee, 2002).Land-cover types were then classified with the estimated harmonic components using an unsupervised classification approach based on a hierarchical clustering algorithm (Lee, 2002).

In this study, the NDVI images were computed for MODIS and GOCI imagery over the Korean Peninsula for 5 years from 2012 to 2016. Section II briefly describes data reconstruction, multi-periodic harmonic modeling and classification algorithm based on a hierarchical clustering algorithm. Section III containsthe results of applying the proposed algorithm to MODIS and GOCI NDVI time series. Section IV has some conclusions.

2. Classification Algorithm

1) Adaptive Data Reconstruction

Satellite data inevitably contain disturbances caused by atmospheric effects and surface anisotropy scattering, which result in unobserved or bad-quality values.They impede estimation of an accurate harmonic model for a certain period of time. In this study, for high-quality data stream, we perform dynamic compositing, which is a feedback system consisting of dynamic synthesis, Bayesian MAP estimation and adaptive harmonic model estimation (Lee, 2002). Dynamic synthesis is a process of continuously recovering unobserved or badquality data using an adaptive quadratic polynomial model. The predicted value of the current time is estimated by using the adaptive polynomial model which showsthe trend ofthe data until the present time, and the bad-quality data is restored by the fitting. The reconstructed data in the dynamic synthesis is used to estimate the adaptive harmonic model after spatial correlation is removed by Bayesian MAP estimation (Lee, 2008). The Bayesian MAP estimation can be interpreted as a process of removing the spatial correlation ofthe observed data to performthe estimation process ofthe time series model independently for each pixel, and it can be said to be the MAP estimation process of the original value using the spatial context information. The dynamic synthesis is a process of restoring the observation values of unobserved or badquality data using short-term time series information. Bayesian MAP estimation is a process of estimating original values considering spatialfactors, and adaptive harmonic model estimation is a process of correcting the Bayesian MAP estimate value by temporal factor.

The adaptive reconstruction system can generate a completely reconstructed image in real time at a given unit period. In addition, it is possible to analyze the characteristics ofsurface vegetation change in real time using the harmonic factor estimated every unit period.

2) Multi-periodic Harmonic Model

Land surface shows complicated physical processes, which are difficult to describe with a simple sinusoidal model (Lee, 2008). A complex model of multiple periods would be more properto represent inter-annual and inner-annual variations of surface parameters, resulting in characterizing a complex curve.Thisstudy proposes a multi-periodic harmonic model, which is expressed as the sum of a series of sine waves.

An image is considered as a set of n pixels and the intensity process can be represented in the following form,

\(Y_t = X_i + \varepsilon_t \)        (1)

\(X_t = m_t + \delta _t ={\mu + \sum {_{k=1}^p}R_{ik}\sin (\omega _{ik}t + \phi_{ik}),i}\in I_n\)        (2)

where I_n = {1, 2, Λ, n} is a set of pixel indices, Y_t and X_t are the observed and original intensity vectors respectively, εt isthe spatially-correlated random noise vector, mt mt is the mean intensity vector associated with local texture, δ_t is the deviation vector from the mean intensity, and p is the number of harmonic periods.

The sinusoid form of Eq. (2) for each pixel is individually restated without the pixel index as

\(x_t =\mu +\sum {_{k=1}^p}R_k \sin(\omega _k t + \phi _k) = \mu +\sum {_{k=1}^p} \alpha _k \cos \omega _k t + \beta _k \sin \omega _k t\)        (3)

Where \(\sqrt {\alpha ^2 _k + \beta ^2_k} \qquad\\ \phi =\tan ^{-1}{\alpha _k \over \beta _k} \ for\ \beta _k \neq 0,\) 

For β_k = 0, øk has a value of 0, π, -π (Bloomfield, 1976). The mean level μ, frequency ω_k, amplitude R_k and phase ø_k are the harmonic components. The parameters of the harmonic model are derived from the temporal trajectory of each pixel’sintensity and can be adaptively estimated by using the exponentially weighted least squares criterion (Lee, 2002). The classification is based on multi-temporal features over the observed area using the estimates of the harmonic components.

3) Multistage Hierarchical Clustering Classification

Vegetation on the earth’s surface varies with time. The time series of vegetationwith seasonal characteristics can be modeled by using four factors of harmonic function: mean, amplitude, frequency and phase. In terms of the temporal characteristics in the vegetation process, the average values of the overall greenness are the mean intensity level over the period, the amplitude reflects the range (size of the wave), frequency is the seasonal periodicity, and the phase represents the initiation time of variation. Vegetation types can be distinguished by the parameters of the harmonic components. Thus the estimated harmonic components are used as feature vectors for land cover classification.

In this study, a multi-stage hierarchical clustering techniquewas employed for unsupervised classification. Hierarchical clustering is one of the most appropriate approaches for unsupervised analysis, which stepby-step merges small clusters into larger ones using similarity (or dissimilarity) coefficient (Anderberg, 1973).In image classification techniques, it is necessary to consider an essential structural characteristic of the information hierarchy (Tanimoto and Klinger, 1990). In the hierarchicalstructure, more than one sub-regions in the lower levels can be merged into large homogeneous regions in the higher levels and this processisrepeated atsuccessively higherlevels. Under the assumption of the hierarchical structure, it is possible to determine natural image segments by combining hierarchical clustering (Lee, 2001).

The multi-stage algorithm consists of two stages. In the first stage, local segmentor, region-growing segmentation is performed by employing a hierarchical clustering procedure with the restriction that pixels in a cluster must be spatially contiguous. In the second stage, global segmentor, hierarchical clustering is carried outforthe segmentsresulting from the previous stage withoutspatial constraintsfor merging.The local segmentation can be considered to be a stage to reduce the obscurity in the image pattern, whereas the global segmentation a classification stage in which the image is grouped into meaningful regions.

Suppose that the image is partitioned in m regions in h level of the multistage hierarchical clustering. Let J_m = {1, 2, …, m} be a set of region indices associated with a G^h = {G_j^h , j ∈ J_M} partition and be a partition at the hth step where G_j^h is a set of the pixels pertaining to region j. For convenience, the index h is omitted and the variable or sets related to the merge of G_r and G_s are indexed by u. That is,

\(G_r \cup G_s = G_u , r,s\in J_m\)        (4)

G_u means a partition state where G_r and G_r in G are replaced by G_u.

Let Z = {Z_i, i ∈ I_n} and Z_i be a feature vector that represents the ith pixel. The similarity coefficient is then obtained in the following:

\(\lambda _u = \triangle _q (Z|G) - \triangle _q (Z|G_u)\)        (5)

where P (Z | G) ∝ exp {–Δ_q (Z | G)} and P(G_u | Z, G) ∝ exp {λ_u}. Therefore, the clustering approach using the similarity coefficient of Eq. (5), λ_u, yields the clustering state with the maximum likelihood condition among all the possible regional configurations at every level. The Δ_q term can be generally represented by a function of the difference between the mean value and the observed value of the feature vector (Lee, 2002).

3. Applications

In this study, the proposed algorithm was applied to MODIS and GOCI NDVI sequential images of daily intervals over the Korean Peninsula for five years from 2012 to 2016. In MODIS, the actual area of the pixel is 250 m × 250 m and one image has 2590 × 4380 pixels. A single GOCI image has 1210 × 2190 pixels, each having a size of 500 m × 500 m.

First, the observed time series data wasrecovered by the adaptive reconstruction system into high-quality time series data stream. The adaptive harmonic model uses the adaptive index to forget the past trends and use only information that is close to the present. In the dynamic synthesis, we use a composite index that indicates the reliability of the observed values. The adaptation or composite index has a value between 0 and 1, and as the value is getting close to 1, the estimation is based on the past trend than the present values. Since the dynamic synthesis is aimed at reconstruction of observed values, a relatively small composite index is adopted to only use information for a short period of time that is closer to the present than a trend over a long period of time. As the adaptive harmonic model is used to estimate the original value, a relatively large adaptation index is employed to use the information over a long period oftime.The adaptive index uses a constant value at every point of time, but the composite index is determined differently by estimating the reliability of the observed value at each point in time. Fig. 1 visually showsthe data reconstruction effect. It is noticed that the temporal profile is represented more clearly after data reconstruction and it has better harmonic model fitting compared with wavelet-based filtering which is one of the represen tative noise reduction methods in frequency domain

Fig. 1. Data reconstruction: blue- original MODIS NDVI data, red- harmonic model fitting after data reconstruction by the proposed method, black- harmonic model fitting after data reconstruction by a wavelet-based method.

Here, for adaptive harmonic analysis of the Korean peninsula, four adaptive harmonic models were used: a one-year cycle (period), a six-month cycle (period), a three-month cycle (period), and a one-month cycle (period). The parameters of harmonic components are estimated from the reconstructed NDVI time series at the end of each year. The mean, minimum and maximum values of amplitude estimated for the reconstruction series of MODIS and GOCI NDVI are shown in Table 1 and 2. As shown in the results of Tables, the one-year period has bigger influence than the factors of other periodssince the Korean peninsula is a distinct four-season region.

Table 1. Estimated amplitude of harmonic model with 4 Periods for MODIS data of 5 years

Table 2. Estimated amplitude of harmonic model with 4 Period for GOCI data of 5 years

We selected meaningful elements to classify landcover types among the estimators of multi-periodic harmonic components as a feature for classification.In thisstudy, a feature vectorforthe i-th pixel is composed as follows,

\(f_i = \begin{bmatrix}\mu \\ R_1 \\ \phi_1 \\ R_2\end{bmatrix}\)

where μ is the mean value, R₁ is the amplitude of the one-year period, ø₁ isthe phase of the one-year period, and R₂ is the amplitude of the half-year period.

The classification results of MODIS NDVI data for 5 years from 2012 to 2016 are shown in Fig. 2(a)-(e). Classification map is consisted ofseven classes and the percentage shown in the legend represents the area occupied by each class. Fig. 2 also showsthe harmonic patterns of vegetation for each class. Table 3-1 to 3-5 contains the estimated values of the harmonic parameters of each class. Classifications results show similar distributions of vegetation each year. In this classification, the first class (Class 1) corresponds to the urban area, and the northern part of the peninsula hasthe second, third and fourth classes(Class 2, 3, and 4). The vegetation types of the fifth, sixth and seventh classes(Class 5, 6, and 7) are mainly distributed in the southern part of the Korean peninsula. The average intensity level in the southern region is higher than in the northern region, but the northern region has higher amplitude. Similarly, the classification results of GOCI NDVI data for 5 yearsfrom 2012 to 2016 are shown in Fig. 3(a)-(e). Table 4-1 to 4-5 contains the estimated values of the harmonic parameters of each class. Comparing these results with those of MODIS, the distribution of vegetation types doesn’t show any significant difference, but the mean intensity level is lower in GOCI and the amplitude is slightly higher in MODIS.

Fig. 2. Classification Map using 4 Harmonic Components for MODIS of Year 2016-2012.

Fig. 3. Classification Map using 4 Harmonic Components for GOCI of Year 2016-2012.

Table 3-1. Estimated Class Harmonic Parameters for MODIS of Year 2016

Table 3-2. Estimated Class Harmonic Parameters for MODIS of Year 2015

Table 3-3. Estimated Class Harmonic Parameters for MODIS of Year 2014

Table 3-4. Estimated Class Harmonic Parameters for MODIS of Year 2013

Table 3-5. Estimated Class Harmonic Parameters for MODIS of Year 2012

Table 4-1. Estimated Class Harmonic Parameters for GOCI of Year 2016

Table 4-2. Estimated Class Harmonic Parameters for GOCI of Year 2015

Table 4-3. Estimated Class Harmonic Parameters for GOCI of Year 2014

Table 4-4. Estimated Class Harmonic Parameters for GOCI of Year 2013

Table 4-5. Estimated Class Harmonic Parameters for GOCI of Year 2012

Four parameters characterize the seasonal variation in temporal profile in a harmonic model.Therefore, we can classify the land cover by the parameter of interest to understand its related characteristic. For instance, Fig. 4 shows mean intensity level and its corresponding classification map for MODIS and GOCI NDVI data in year 2016, which shows the distribution of land-cover types in the overall greenness’s point of view. Similarly, amplitude of one-year period and its corresponding classification map, and phase of oneyear period and its corresponding classification map, respectively, are displayed. This approach based on each parameter makesit possible to analyze the overall vegetation volume according to the average intensity level, the peak value at the peak of growth through amplitude, and the starting point of the growth phase through phase during the analysis period.The selection of number of classis depending on which class can be classified with the parameters of interest.

Fig. 4. Examples of estimated parameters and their corresponding classification maps for MODIS and GOCI NDVI of year 2016.

Change detection based on classification can be also produced in the proposed approach. For instance, Fig. 5 and 6 show the change detection maps according to the total amount of greenness of MODIS and GOCI NDVI data during 2012 and 2016, respectively, which isshown as annual average level difference (increase / decrease). The black area isthe area with little change, the red area is the area where the incremental change occurred, and the blue area is the area where the decreased change occurred. Similarly, the change analysis can be applied to each component of the harmonic model such as amplitude and phase. In the case of amplitude, the change detection map can be made with change types, the area with little change (difference (-π) / 60 and π / 60), the area where the incremental change occurred (difference is greater than π / 60), and the area where the decreased change occurred (the difference is less than -π / 60) and so on. This approach can be performed for monitoring vegetation changes per specific parameter such as phase and specific period such as annual or semiannual. In this way, it is possible to understand the change ofland cover according to the characteristics of the vegetation related to the harmonic parameters.

Fig. 5. Annual difference in mean level in MODIS NDVI data during 2012 and 2016.

Fig. 6. Annual difference in meanz level in GOCI NDVI data during 2012 and 2016.

4. Conclusions

Classification is a very useful technique for extracting land-cover information in analysis of remote sensing data. Various classification methods and approaches have been developed. An appropriate classification method can be selected depending on the applications. This study focused on developing a multi-temporal classification method for a sequential ofremote sensing data collected regularly at short time interval. Since multi-temporal data include rich temporal dynamic information over time, it is appropriate to analyze the data in time domain using featuresthat can explain the temporal process.This approach is an alternative to the conventional techniques solving the problem of signal variability that occurs in spectral analysis.

In this study, a classification method of analyzing temporal variability based on the harmonic model is proposed including data reconstruction as a preprocessing for high quality data stream.The experimental results say that harmonic-model based approach may be the most efficient and plausible temporal technique to analyze time series data exhibiting seasonal trends such as MODIS and GOCI NDVI representing landcover vegetation activity. Land-cover classification can be performed with a feature of interesting characteristics.This enables usto understand land cover and land-cover changes in different point of view.

Fortunately,mostofvegetationprocessesshows annual behavior and inter-seasonal phonological variability. Each vegetation type can be characterized by a unique temporal profile to the extent of its development cycle. Therefore, harmonic-model based classification reflects different sources of temporal variation by using the estimated values of the parameters associated with the temporal profile observed in the sequential images. Classification results can also be utilized for change detection. The change detection of the earth’s surface is to be done timely and accurately and then the relationship and interaction between natural phenomena and humans can be better analyzed and understood as a result. The proposed method can be also utilized in change detection (Jung and Jang, 2014).