SPECTRAL UNMIXING ANALYSIS OF TIME SERIES LANDSAT 8 IMAGES

Temporal analysis of Landsat 8 images opens up new opportunities in the unmixing procedure. Although spectral analysis of time series Landsat imagery has its own advantage, it has rarely been studied. Nevertheless, using the temporal information can provide improved unmixing performance when compared to independent image analyses. Moreover, different land cover types may demonstrate different temporal patterns, which can aid the discrimination of different natures. Therefore, this letter presents time series K-P-Means, a new solution to the problem of unmixing time series Landsat imagery. The proposed approach is to obtain the “purified” pixels in order to achieve optimal unmixing performance. The vertex component analysis (VCA) is used to extract endmembers for endmember initialization. First, nonnegative least square (NNLS) is used to estimate abundance maps by using the endmember. Then, the estimated endmember is the mean value of “purified” pixels, which is the residual of the mixed pixel after excluding the contribution of all nondominant endmembers. Assembling two main steps (abundance estimation and endmember update) into the iterative optimization framework generates the complete algorithm. Experiments using both simulated and real Landsat 8 images show that the proposed “joint unmixing” approach provides more accurate endmember and abundance estimation results compared with “separate unmixing” approach.


INTRODUCTION
Hyperspectral remote sensing technology consists of acquiring a set of images capturing a spatial scene at a few hundreds of contiguous spectral bands.However, the significant spectral information they convey is somewhat compromised by their lower spatial resolution.This limitation, combined with the complexity of the ground targets and environmental conditions, lead to the observed pixel spectra composed of several pure materialsreferred to as endmembers.Spectral unmixing aims at identifying a number of endmembers and their abundance fractions in each pixel.
So far, spectral mixture analysis has been widely used to characterize the spectral signatures of land cover types at a single time point.The limited spectral information possessed by each time frame of the image is not adequate to accurately estimate endmembers.Consequently, spectral unmixing still presents a challenge, as different ground components may demonstrate similar spectral signatures in a remote sensing image due to the low spectral resolution of multispectral imagery.Therefore, multitemporal analysis of remote sensing images has drawn increasing interest in recent years, mainly for change detection, feature selection and classification.For instance, Du et al. (Du et al., 2005) described a change detection approach based on a linear mixing model for multi-temporal CASI (Compact Airborne Spectrographic Imager) images over crop fields.To date, temporal unmixing analysis has been mostly conducted using the time series of MODIS or Hyperspectral images (Diao and Wang, 2016).Given above consideration, spectral unmixing a sequence of Landsat 8 images captured over the same area can be a significant interest.
Several attributes of Landsat 8 are wide scope of coverage, higher spatio-temporal resolution and cost-free status.More importantly, these data are well available at regular time intervals.But the unmixing analysis of a time series of Landsat imagery has rarely been studied.Dedicated methods in the field of multitemporal spectral unmixing have only begun to emerge (Iordache et al., 2017).Somers et al. (Liu et al., 2012) proposed a method to minimize the effect of endmember variability in the image by grouping multiple spectral signatures (or bundles) to describe a particular endmember class in the scene.Although all these methods have their own respective advantages, the acquisition of the temporal signatures of all the endmembers in the images is still a challenge.Spectral unmixing would be more straightforward if we have "pure" pixels.Therefore, we intend to obtain the "purified" pixels from the mixed pixels in order to achieve a favorable unmixing performance.The definition of a "purified" pixel is that the residual of the mixed pixel after excluding the contribution of all nondominant endmembers.Accounting for the temporal information, we will refer to as time series K-P-Means algorithm, which alternates iteratively between two main steps (abundance estimation and endmember update) until convergence to yield final endmember estimates (Xu et al., 2014).Endmember update refers to utilize the endmembers that have already been updated in the current iteration.The superiority of the analysis of time series of unmixing with the K-P-Means algorithm is evaluated by experiments on both simulated and real multispectral images.
The paper is organized as follows.The proposed time series K-P-Means model accounting for temporal information is introduced in Section 2. Experimental results obtained on synthetic and real data are reported in Section 3. Section 4 finally concludes this work.

Problem Formulation
A linear mixing model (LMM) is adopted in this letter since it is appropriate to describe remote sensing data when the declivity of the scene and microscopic interactions between the observed materials are negligible (Bioucas-Dias et al., 2012).The LMM consists of representing each image pixel by a linear combination of the endmember matrix and abundance matrix with the Gaussian noise.The model can thus be written as where xi is a P × 1 dimensional image pixel, s ik is a 1 × 1 nonnegative abundance vector, and a k denotes a P × 1 matrix.Finally, n represents an additive noise.In order to achieve an efficient endmember estimation, it is reasonable to take advantage of the "purified" pixels that are only due to dominant endmember.We refer to xi after removing the contribution of nondominant endmembers as "purified" pixels.
It is desirable to utilize the good abundance information to obtain "purified" pixels for elegant endmember extraction.On the other hand, the accurate endmember estimation also can enhance the accuracy of abundance estimation.Therefore, spectral unmixing is interpreted as an optimization issue, which is solved by iteratively alternating the estimation between the endmember update and abundance estimation.The following section describes the time series K-P-Means model.

time series K-P-Means model
In this section, we consider Landsat 8 images acquired at T different time frames over the same scene, assuming that K endmembers are present in the resulting time series.Given an a priori known number of endmembers K, the time series K-P-Means model is based on K-P-Means algorithm, which characterizes a class by the dominant endmember, whose fractional abundance is the biggest (Xu et al., 2014).At each time t, the model of time series K-P-Means can be expressed as where t = 1, 2, ..., T and i = 1, 2, ..., N , x k it denotes the ith image pixel at time t, art is the rth endmember at time t, sir is the proportion of the rth endmember in the ith pixel.Finally, n is independent and identically distributed Gaussian noise.The proposed time series K-P-Means accounting for temporal information is written as T is a (P ×T )×1 matrix containing the ith pixel of the all time frames, Ar = [(ar1) T , (ar2) T , ..., (arT ) T ] T denotes a (P × T ) × 1 matrix containing the rth endmember of the all time frames, sir denotes a 1 × 1 abundance.The object function of time series K-P-Means can be formulated as where Yi is the "purified" pixel that is only due to the contribution of dominant endmember, therefore Yi can be formulated as Following ( 4), where l are the labels of pixels, given {Ar}, pixel labeling requires solving the following optimization issue: The above equation means that Xi is associated with the kth endmember A k , which will take the largest coefficient s ik when the representation error is minimized (Xu et al., 2014).

Abundance estimation
According to the previously described model, one essential step is to estimate {sir} given {Ar}, which is a nonnegative least square (NNLS) issue, and is typically formulated as (Bro and De Jong, 2015) (Lawson and Hanson, 1995): Based on this above formulation, a popular approach for NNL-S is an active-set method, which was proposed by Lawson and Hanson in (Lawson and Hanson, 1995) and modified by Bro and Jong for fast computation (Bro and De Jong, 2015).

Endmember estimation
Another step in K-P-Means algorithm is the estimation of , with the following equation: The complete algorithm will be tested on the simulated and real multispectral images.

Complete algorithm
Assembling two main steps (abundance estimation and endmember update) into the iterative optimization framework generates the complete algorithm of time series K-P-Means.The entries of input are the spectral stack X, number of clusters K and a predefined maximum number of iteration.First, the vertex component analysis (VCA) is used to extract endmembers for endmember initialization.Then, nonnegative least square (NNLS) is used to estimate abundance maps by using the endmember.Finally, the estimated endmember is the mean value of purified pixels.

EXPERIMENTS
We evaluate the performance of time series K-P-Means using the Xi and Ar, hereafter called "joint unmixing".We compare it to the results achieved by using xit and art, which is referred to as "separate unmixing".The results of the experiments are evaluated by spectral angle distance (SAD), spectral information divergence(SID), abundance angle distance (AAD), abundance information divergence (AID).

Experiment on simulated images
The proposed method has been applied to a synthetic time series of ten images, composed of four endmembers of size 58×58 with 100 bands.First, we design an abundance map composed of four real spectra, which are randomly selected from the USGS Digital Spectral Library, sampled on 10 spectral reflectance bands.Then, we set a spectral scale factor as one period of a sinusoid over the ten time instants to simulate the time series of endmembers (Henrot et al., 2016).The mixtures correspond to linear mixture of 4 endmembers.The resulting image has finally been corrupted by adding Gaussian noise.In the simulated study, in order to investigate the algorithm robustness to noise corruption, they are tested on simulated images with different noise levels measured by the signal-to-noise ratio (SNR) (Miao and Qi, 2007).The performances of these two methods are compared at different noise levels in terms of AAD, AID, SAD and SID.In these four statistics, smaller value means better result, as reported in Figure 1.
Overall, the "joint unmixing" approach achieves a better performance than the "separate unmixing" approach across all noise levels, indicating that spectral unmixing analysis on time series Landsat 8 images can achieve a good estimation result and that K-P-Means algorithm is capable of accurately estimating endmembers in highly mixed and noisy circumstance.Note that here, compared to the separate unmixing method, the endmember estimation of joint method measured by SAD and SID seem to be robust with SNR decreasing from 45 to 20.As we can see, the numerical variation of SAD and SID is very low.When SNR =10, we noticed that "joint unmixing" achieves much lower SAD and SID values than "separate unmixing".In terms of abundance estimation, the results display the similar pattern, joint method produces better performance than separate method, although in the high noise circumstance, there was a small discrepancy in AAD and AID.In general, it can be seen from the Figure 1 that the "separate umixing" performs worse than the "joint umixing".

Test on the real images
The time series of Landsat 8 images, which were captured over the scene of the Bohai Gulf, is used to test the proposed algorithms.There are 11 images with minimal cloud cover were collected by OLI and TIRS in 2014 and these are available in the Geospatial Data Cloud website.The temporal analysis was done using the images collected in March 5 and 21, April 6 and 22, May 8, June 9, July 11 and 27, Aug. 12, Sep. 13 and Oct. 15 of 2014.Each image comprises 456 × 453 spatial pixels belonging to 4 different land cover types.Each image has spatial resolution of 30 m and contains 7 wavelengths.
In this section, we compare the separate and joint unmixing approach on this data set.The spectral unmixing results of the "separate unmixing" and the "joint unmixing" on the different time frames are displayed in the Figure 2 and Figure 3.As we can see, the unmixing performances of these two methods seem to show a similar pattern.Especially the "endmember 1" and "endmember 2" seem correctly extracted in both cases.But, combined with the abundance maps, the joint method unmixes "endmember 3" and "endmember 4" better than the separate method.Other time frame results also demonstrate that joint method yields visually better performance.Moreover, accounting for the temporal information in the joint unmixing approach provides a good estimation result, as shown in Figure 4.
In order to quantitatively evaluate the performance of the two methods, we run the comparison for seven different time frames and the value of SAD and SID of the estimated spectra reveal that "joint unmixing" outperforms "separate unmixing", as shown in Table I 1.When SNR=40, performance of the "joint unmixing" approach and "separate unmixing" approach, measured by SAD and SID, over different time frames.In these statistics, smaller value means better result.

CONCLUSION
This paper has presented a solution to the problem of unmixing a time series of Landsat 8 images based on the K-P-Means algorithm.The time series K-P-Means algorithm was to jointly process a time series of multispectral images for an optimal unmixing result.Indeed, sequential endmember estimation from a set of Landsat images captured over the same area could provide improved performance when compared to independent image analyses.Moreover, different land cover types may exhibit different temporal patterns, which could aid the discrimination of different natures.Therefore, temporal analysis could be of significant interest.The purpose of time series K-P-Means algorithm was to obtain the "purified" pixels for enhanced endmember Figure 2. Endmembers obtained by the "separate unmixing" method and "joint unmixing" method.Estimated spectra from the fifth, seventh and tenth time frame.The "separate unmixing" results are displayed in (a), (b) and (c).The "joint unmixing" results are displayed in (d), (e) and (f).Note that the green indicates the true endmembers, the black represents the estimated endmember.The estimated endmembers are closer to the true endmembers.See e.g. the fifth, seventh and tenth time frames of endmember 3 and 4. Figure 3. Abundance maps obtained by the "separate unmixing" method and "joint umixing" method.The "separate unmixing" results are displayed in the first row.The "joint unmixing" results are displayed in the second row.As we can see, the unmixing performances of the two methods on "endmember 1" and "endmember 2" show similar patterns.However, combined with the estimated spectra in the Figure 3, the joint method unmixes "endmember 3" and "endmember 4" better than the separate method. .Endmembers obtained by the "separate unmixing" method and "joint unmixing" method.The "joint unmixing" approach accounts for the temporal information.Note that the green indicates the true endmember, the black represents he estimated endmember.The estimated spectra by joint unmixing are closer to the true endmembers, especially "endmember 4" extracted correctly by the joint method.Other endmember estimation results also indicate that joint method yields better performance.estimation, which could, in turn, improve the accuracy of abundance estimation.Consequently, time series K-P-Means algorithm could be treated as an iterative optimization problem by alternating the estimation between two main steps (endmember estimation and abundance estimation).The performance of the presented approach proved that "joint unmixing" approach outperformed "separate unmixing" approach and that time series K-P-Means was capable of accurately estimating both endmembers and the corresponding abundance.

Figure 1 .
Figure 1.The comparison of "joint unmixing" and "separate unmixing" at different noise levels in terms of (a) SAD, (b) SID, (c) AAD, and (d) AID.In these four statistics, smaller value means better result.
Figure4.Endmembers obtained by the "separate unmixing" method and "joint unmixing" method.The "joint unmixing" approach accounts for the temporal information.Note that the green indicates the true endmember, the black represents he estimated endmember.The estimated spectra by joint unmixing are closer to the true endmembers, especially "endmember 4" extracted correctly by the joint method.Other endmember estimation results also indicate that joint method yields better performance. .