- Email: [email protected]
Using NDVI time series to diagnose vegetation recovery after major earthquake based on dynamic time warping and lower bound distance
Ecological Indicators 94 (2018) 52–61
Contents lists available at ScienceDirect
Ecological Indicators journal homepage: www.elsevier.com/locate/ecolind
Using NDVI time series to diagnose vegetation recovery after major earthquake based on dynamic time warping and lower bound distance ⁎
Xuelei Zhanga, Ming Wanga,b, , Kai Liua,b, Jun Xieb, Hong Xua,
T
⁎
a State Key Laboratory of Earth Surface Processes and Resource Ecology/Academy of Disaster Reduction and Emergency Management, Faculty of Geographical Science, Beijing Normal University, Beijing 100875, China b Key Laboratory of Environmental Change and Natural Disasters, Ministry of Education, Beijing Normal University, Beijing 100875, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Earthquake Vegetation recovery NDVI Dynamic time warping Lower bound distance Difference measurement index
Major earthquake occurred in mountainous areas usually cause large number of landslides that lead to severe impact to local vegetation cover and growth. The negative influence to vegetation may last for many year and vegetation recovery may experience dynamic fluctuation. Existing methods for vegetation recovery diagnosis face difficulty in capturing the dynamic behaviours both within and between years that makes the interannual comparison impossible. This paper proposes a new method to diagnose regional vegetation recovery after a major earthquake by defining a difference measurement index (DMI) using MODIS NDVI time series at 8-day interval. This differs from many existing methods in its quantification of the difference between the studied time series and historical samples, by using a proposed algorithm consisting of lower bound distance and dynamic time warping. This algorithm can better differentiate vegetation disturbance from its natural fluctuation. Second, the method investigates relatively regional vegetation recovery via a dynamic index, the DMI. Vegetation conditions in different years can be compared with a historical benchmark and measured by DMI. This makes it possible to diagnose dynamic vegetation recovery and generate a series of interannual spatial distributions of regional vegetation state.
1. Introduction When major earthquakes strike in mountainous regions, a large number of landslides, debris flows and quake lakes often occur and destroy vast areas of vegetation cover (Liu et al., 2008). After the Wenchuan Earthquake, several studies of regional ecological change and vegetation damage assessments were carried out using satellite data at multi-temporal scales (Ge et al., 2009). Damage to vegetation can be reflected in satellite images, for example via decreased Normalized Difference Vegetation Index (NDVI) values (Wang et al., 2014b). Much research has been undertaken concerning the immediate impact of earthquake-induced landslides on vegetation (Mantovani et al., 1996, Metternicht et al., 2005, Chen et al., 2009, Chong et al., 2009, Cui et al., 2012, Wang et al., 2014b). However, existing methods have difficulty in assessing the dynamics of vegetation recovery. More efforts should be paid to develop new methods for revealing the dynamic process of vegetation recovery and conduct case studies for crosscomparison of various methods. Remote sensing data can effectively diagnose vegetation recovery of mountainous regions (Tømmervik et al., 2003, Turner et al., 2007, Ni et al., 2009, Wang et al., 2009, Lu
et al., 2012), although researchers always face challenges from unsatisfactory data quality (e.g., cloud cover) and natural variability (e.g., changes in phenology) (Lu et al., 2012). The vegetation recovery rate (VRR) (Lin et al., 2006) method compares vegetation changes at a few cross-sectional data before and after a major disturbance. It provides an indicator of vegetation changes over time, however, it cannot give a complete picture of vegetation change since the occurrence of disturbance. The COI (Count of Outlier Index) (Wang et al., 2014b) determines an outlier if the NDVI value is double standard deviation less than past mean at a certain time window. The total number of outliers (the COI value) reflects the overall vegetation recovery status since the occurrence of the major disturbance, but it may generate false detection of outliers due to interannual phenological changes. Another approach based on the time series of satellite data using regression algorithms is to assess the trend of vegetation change and recovery rate (Röder et al., 2008, Van Leeuwen, 2008, Gouveia et al., 2010). However, the change in trends may significantly affected by the insufficient data length and the presence of abnormal values in original data. Very often, filters are used to eliminate the effect of cloud and abnormal points (e.g., cloud cover)
⁎ Corresponding authors at: State Key Laboratory of Earth Surface Processes and Resource Ecology/Academy of Disaster Reduction and Emergency Management, Faculty of Geographical Science, Beijing Normal University, Beijing 100875, China (M. Wang for general correspondence and H. Xu for data). E-mail addresses: [email protected] (M. Wang), [email protected] (H. Xu).
https://doi.org/10.1016/j.ecolind.2018.06.026 Received 23 October 2017; Received in revised form 6 June 2018; Accepted 12 June 2018 1470-160X/ © 2018 Elsevier Ltd. All rights reserved.
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 1. Geographic location of study area (eastern Pingwu County).
(Cuevas-González et al., 2009). Savitzky–Golay filter (SG filter) is widely used based on a weighted moving average algorithm (Chen et al., 2004), but the determination of filter parameters usually influences analysis results (Chen et al., 2004, Diouf et al., 2015, Yang and Qi, 2017). Taking the maximum values of every three consecutive bands (Max3 cloud filter) is another effective filter (Yang et al., 2013), but frequent and heavy cloud may not be eliminated effectively. In processing time series data, we need to eliminate obvious abnormal signals (e.g., cloud cover) while retaining signals of sudden change that may be induced by vegetation damage. Therefore, careful processing is required. In this paper, we develop a new method based on 8-day interval MODIS (Moderate Resolution Imaging Spectroradiometer) NDVI time series to diagnose vegetation recovery after a major earthquake. This differs from many existing methods in the way of quantifying the difference between a time series and historical samples by using a proposed DTW-LBD algorithm that consists of LBD (lower bound distance) and DTW (dynamic time warping). LBD is used for increasing the substantial difference between time series while neglecting small distinctions (i.e., NDVI variation caused by moderately increased or reduced seasonal rainfall), and time warping (TW) similarity measure that allows sequences to be stretched or compressed along the time axis. The DTW algorithm (Keogh and Ratanamahatana, 2005) allows elastic shifting of the time axis to accommodate sequences that are similar but out of phase. DTW is used for eliminating any errors induced by interannual phenological change (i.e., growing season delay due to extended winter period). The DTW-LBD algorithm therefore can better differentiate vegetation disturbance from its natural fluctuation and enables diagnosis of dynamic vegetation recovery and generation of a series of interannual spatial distributions of regional vegetation state. The objectives of this paper are to, first, develop a new method that can assess the dynamic change of vegetation recovery after major earthquake by minimizing the influence of interannual phenological changes and persistent small difference, and second, use a case study of the Wenchuan earthquake stricken area to verify the validity and applicability of the proposed method.
2. Materials and methods 2.1. Study area The Wenchuan region in southwest China is ecologically diverse and but fragile. At 2:28p.m. local time on May 12th, 2008, a magnitude 8.0 earthquake occurred in Sichuan Province, followed by thousands of associated aftershocks, resulting in more than 200,000 landslides over an area of 110,000 km2 (Wang et al., 2014b). Eastern Pingwu County was one of the most severely affected areas during the Wenchuan Earthquake, with Modified Mercalli Intensity scale IX and X. Southeast monsoon dominants the local climate. Most rainfall concentrate from July to October with annual average precipitation of 800 mm and a standard deviation of 167 mm (Wang et al., 2014a, Li et al., 2018). The average temperature is 14.8 °C with daily minimum of −4 °C and maximum of 39.1 °C. This geographic extent shown in Fig. 1 was selected to test the proposed DTW-LBD algorithm for vegetation recovery diagnosis. 2.2. Materials The MODIS sensor is on board the Terra and Aqua satellites. We used MODIS MOD09Q1 products (Wang et al., 2014b), with spatial resolution 250 meters and temporal resolution nearly 8 days. MOD09Q1 contains both red light (red) and near-infrared (NIR) bands. We converted the projection to WGS-84, UTM48N according to the location of the study area (Wang et al., 2014b). In that area, 690 images were acquired from January 1, 2001 to December 31, 2015. The NDVIs for all cells were calculated based on reflectance in the NIR and red bands (Zhang et al., 2010, Wang et al., 2014b)
NDVI = (NIR−red )/(NIR + red )
(1)
Cloud, aerosol, and vapor in the atmosphere usually cause an abnormal decrease of the NDVI, thereby not reflecting the actual vegetation condition. Residues that were not detected during the masking process of atmospheric noise also generate abnormal valleys in the time
53
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
series. Here we use Max3 cloud filter with outliers deletion to achieve this goal by taking the following steps. Firstly, per reference (Yang et al., 2013), we took the maximum value of every three consecutive bands to form a new time series (Max3 cloud filter):
NDVIi = max (NDVIi − 1, NDVIi , NDVIi + 1, )
(2)
In this formula, i means the ith time band in original NDVI time series. Secondly, we remove the NDVI values higher than 0.99, and use the average of the left and right NDVI values to substitute the abnormal value. Thirdly, we define an abnormal valley when the NDVI experiences a sudden reduction followed by an immediately sudden recovery. If the value of the previous NDVI is 0.2 higher than the current point while the value of the next NDVI is also 0.2 higher, we define this point as an abnormal valley. In addition, there is a phenomenon that abnormal valley may have more than one point. We also use the average of the left and right NDVI of the abnormal valley to substitute the values in the valley.
Fig. 2. Typical cases of LBD calculation.
minimum difference between two sequences, DTW maps each element of a sequence to one or more neighbouring elements of another sequence. Given any two non-null sequences A and B , the time warping LBD is defined as follows.
2.3. Grid representation and LBD The time series data can be divided by predetermined grid cells and transformed into a form of grid representation (Sun, 2012). A = {a1, a2 , ⋯, an} is the original time series and ai is the value at time i . If the grid width (Wg) is given, the grid representation of A can be defined as GRA = {〈AR1 , ARB1 〉, 〈AR2 , ARB2 〉, ⋯, 〈ARn , ARBn 〉} , in which ARi is the row of ith value of A in the grid and ARBi represents which part of a grid into which the ith value is entered. ABi is an integer denoting the rower. If the ith value of A is above the center line of a grid, ARBi = 1; otherwise ARBi = 0 . The centre line of a grid is calculated by the average of upper and lower bound. In contrast with the widely used Euclidean distance, LBD is used to quantify the difference between two time series A and B based on their grid representations GRA = {〈AR1 , ARB1 〉, 〈AR2 , ARB2 〉, ⋯, 〈ARn , ARBn 〉} and GRB = {〈BR1, BRB1 〉, 〈BR2 , BRB2 〉, ⋯, 〈BRn , BRBn 〉} .
DTWLBD (A, B ) =
Here, First (A) denotes the ‘current’ first element in A, and Rest (A) denotes the sequence consisting of the elements in A other than First (A) . This equation shows an iteration procedure in which the ‘min’ term repeats the DTW–LBD algorithm and the First ( ) and Rest ( ) functions keep updating the first and remaining sequences (Zhu and Shasha, 2003).
n
2.5. Measurement of vegetation recovery
n
⎨ ∑ (AR −BR −1)2Wg 2 , ifAR ≠ BR , i i i i ⎪ ⎪ i=1 ⎪ and (ARi −BRi )(ARBi −BRBi ) ⩽ 0 ⎪ ⎩ 0, if ARi = BRi
2 ⎧ DTWLBD (A, Rest (B )) First (B )) + min DTWLBD 2 (Rest (A), B ) ⎨ 2 ⎩ DTWLBD (Rest (A), Rest (B ))
(4)
⎧ 2 2 ⎪ ∑ (ARi −BRi ) Wg , ifARi ≠ BRi , ⎪ i=1 ⎪ ⎪ and (ARi −BRi )(ARBi −BRBi ) > 0 LBD (A, B ) =
LBD 2 (First (A),
Based on the proposed DTW–LBD algorithm, a difference measurement index (DMI) can be defined: ym + x
DMIx (k ) =
(3)
Fig. 2 uses five typical cases to illustrate the rules in formula (3). When A and B is on the same grid as showed (case (a) and (b)), the LBD of A and B is zero. When A and B is on different grids, there are three typical cases in which the difference between A and B is approximately estimated by counting the grids between A and B. Using the LBD algorithm, the difference of the ith observed values is decreased. In case (c) and (d), the LBD is estimated to be one grid (Wg) as (ARBi −BRBi ) ⩽ 0 , while the LBD is two grids (2 Wg) for case (e) as (ARBi −BRBi ) > 0 .
∑y = y1 DTWLBD (Tyi, Tyj ) Cm2 + 1
ym
−
∑y = y1 DTWLBD (Tyi, Tyj ) Cm2
(5)
Here, DMIx (k ) denotes the vegetation recovery status at location k ym in year x . The term ∑y = y1 DTWLBD (Tyi, Tyj )/Cm2 represents the average DTW–LBD of any two time series Tyi and Tyj in samples from year y1 to ym . This is considered a benchmark of DTW–LBD obtained from historical samples that were not influenced by disturbances (e.g., earthquakes). In this article, y1 = 2001 and ym = 2007 . The term ym + x ∑y = y1 DTWLBD (Tyi, Tyj )/Cm2 + 1 denotes the average DTW–LBD of any two time series Tyi and Tyj in samples from year y1 to ym plus year x . In this article, x = 2008, 2009, ⋯, or 2015. This term differs from the benchmark by adding one more time series of the investigated year. Therefore, the difference between these two terms is larger if the added time series from year x is significantly different from the historical average. DMIx (k ) takes the DTW-LBD difference comparing to historical average. This makes DMIx (k ) comparable between years. As a result, the bigger DMI is, the worse vegetation recovered. Linear regression can be done to analyze the difference measurement trend (DMT) based on calculated annual DMIs, as follows.
2.4. DTW of time series with LBD There is always natural interannual fluctuation of vegetation. Time series patterns may change and lead to dramatic differences if calculated using simply the distance at the same time. Here, we use a time warping (TW) similarity measure that allows sequences to be stretched or compressed along the time axis. The DTW algorithm (Keogh and Ratanamahatana, 2005) allows elastic shifting of the time axis to accommodate sequences that are similar but out of phase. To find the
DMIx (k ) = DMT × x + β0 54
(6)
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Step 1: Prepare MODIS images and use the NIR and red bands to calculate the NDVI time series. Step 2: Apply the filtering algorithm to eliminate possible cloud effects. Step 3: Determine a grid size and transform the NDVI time series into grid format, including rows and 0/1, denoting locations on a grid. Step 4: Calculate DMI for all locations based on the historical benchmark samples and samples from the investigated year, using the DTW–LBD algorithms. Step 5: Calculate DMT by using the linear regression based on DMI results. Step 6: Perform DMI and DMT mapping to illustrate vegetation status and recovery trends.
Table 1 High-resolution satellite image in study area. Sensor
Date
Resolution (panchromatic)
IKONOS
November 2, 2002 and November 8, 2002 April 23, 2010 and August 15, 2010 January 29, 2014 August 20, 2015
1m
WV02+qb ZY-3 GF-2
0.5 m 2.5 m 0.8 m
Here, β0 and β1 are the regression coefficients. If DMT < 0 , the DMI at the analysed location decreases with the increase of years and vegetation tends to be more similar to the historical benchmark, implying a positive trend of recovery. If DMT > 0 , the vegetation tends to have a more significant difference relative to the historical benchmark, implying a negative trend of recovery. Therefore, the distribution of regional recovery trends can be obtained by using values of DMT , and it can be used for qualitative assessment of the recovery trend.
3. Results 3.1. Effect of DTW-LBD We chose some typical NVDI time series from the study area to demonstrate the effectiveness of the proposed method and compare with other methods. As shown in Fig. 3, the red line is one of the NDVI time series of 2001, and the blue line is one of the NDVI time series of 2014. The original NDVI time series have big difference due to the influence of cloud, and after application of the cloud filter, the difference is mainly related to the interannual phenological change. In term of vegetation condition, these two NDVI curves in 2001 and 2014 are considered to be similar as no significant difference is observed other than the phenological change. If we use difference methods to calculate the distance between these two NDVI time series, the results can be very different. As shown in Table 2, Euclid distance is highest. The adjusted Euclid distance (Yang and Qi, 2017) (Adjusted ED) applied the SG filter and a cross-correlogram spectral matching algorithm to eliminate the influences of cloud and interannual phenological changes, but still it is much higher than the
2.6. Accuracy assessment For assessing accuracy of this new method, we used some high-resolution satellite images (Table 1) and statistics analysis to verify the assessment results. We also plotted the NVDI distributions and did statistics analysis of their mean and variance of those picked areas in different years. Z-test and F-test are used to check whether significant difference exists in the mean and variance of NDVI among different years, respectively. 2.7. Step-by-step procedure The following steps were used to generate maps of vegetation recovery in various years, based on the MODIS NDVI time series.
Fig. 3. Influence of interannual phenological change.
Table 2 Comparison of distance calculated using different methods.
Case with phenological change (Fig. 3) Case with persistent small difference (Fig. 4)
Euclid Dist.
Adjusted ED
Dist.w/DTW
DTWLBD (Wg = 0.01)
DTWLBD (Wg = 0.05)
DTWLBD (Wg = 0.1)
1.8789 1.5803
1.1012 0.6003
0.4547 0.1784
0.4144 0.1371
0.3122 0.0500
0.2236 0
55
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 4. Influence of persistent small difference.
albeit with some fluctuation. Many non-landslide areas showed a negative recovery trend, indicating that earthquake-stricken areas still faced major challenges because of subsequent geo-hazards or man-made disturbances.
proposed distance with DTW and LBD. If we use NDVI time series in 2001 as the base pre-seismic NDVI time series and the distance between these two NDVI time series to represent vegetation recovery, the Euclid distance and Adjusted ED may not provide an accurate assessment. So DTW could eliminate the influence of interannual phonological changes. Persistent small difference between two time series can also be accumulated into a big difference. As shown in Fig. 4, the red line is one NDVI time series of 2003, and the blue line is one NDVI time series of 2006. There is persistent small difference between these two NDVI time series, and obviously the vegetation condition in these two years is considered to be similar. As shown in Table 2, Euclid distance is highest. With the use of LBD algorithm, this kind of persistent small difference can be eliminated.
3.3. Accuracy assessment The DTW-LBD method is compared with the COI method proposed by Wang et al., (2014b) and the adjusted ED method proposed by Yang and Qi (2017). We picked two areas for the purpose of verification: Ma An Shi landslide (I in Fig. 6) and a typical forest without disturbance (VII in Fig. 6). From the high-resolution satellite images of 2010 in Ma An Shi landslide (Fig. 7), we can see that the landslides caused vegetation deterioration. The vegetation was recovering gradually in the following years such as 2014 and 2015. The mean of NDVI is increasing. In 2015, the mean NDVI is even higher than that in 2002 (significant at p < 0.05), but the variance is also higher (significance at p < 0.05), implying that many local areas experienced vegetation degradation although the overall vegetation was recovering. The results of DMI well demonstrated this dynamic recovery process. The DMI results were compared with the COI and adjusted ED results, suggesting that the DMI results in general agreed with the COI results while the COI results may sometimes exaggerate the difference. The adjusted ED results did not provide sufficient differentiation in terms of spatial patterns. We also selected an area of typical forest that experienced no obvious disturbance (Fig. 8). The mean and variance of the NDVI in 2014 and 2015 remain similar (significance at p < 0.05), which agreed with the DMI results. However, the other two methods exaggerated the difference in some years, which may not accurately reflect the actual situations.
3.2. Vegetation recovery diagnosis We used 690 images in the study area to create the NDVI time series, spanning 2001–2015. By taking samples from 2001 to 2007 as the benchmark, annual DMIs were obtained for 2008–2015. Fig. 5 shows that the study area had dynamic change of vegetation recovery after the major earthquake in 2008. DMIs in all eight figures are at the same scale, so interannual differences are comparable. Large DMI values in 2008 correspond to locations of large landslides, mainly triggered by the earthquake. That means severe vegetation disturbances right after the major earthquake concentrated in the areas with large co-seismic landslides. The heaviest damage to vegetation (in red or orange colour) occurred in the northern part of the area, where two top ten largest landslides are located. These include Ma An Shi (west) and Zheng Jia Shan (east). Starting from 2009, the vegetation in those landslides areas recovered gradually, while vegetation disturbance started to appear in those non-landslide areas, particularly in the southern part of the study area. The vegetation recovery in the study area shows an overall fluctuating and spatial shift. Fig. 6 shows the recovery trend of the study area for the seven years after the Wenchuan earthquake. All large-scale landslides in the study area (I-VI circled in green) had a positive recovery trend, and vegetation of Shi Ban Po landslide, Wen Jia Ba landslide, Zheng Jia Wan landslide recovered (significance at p < 0.1). These landslide areas had heavy vegetation damage in 2008 but gradually recovered their vegetation,
4. Discussion 4.1. Effect of cloud filter selection Savitzky–Golay filter (SG filter) is a widely used method to resolve cloud noise in NDVI time series using a weighted moving average filter with the weighting given as a polynomial (Chen et al., 2004). The key
56
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 5. Assessment of annual vegetation change after Wenchuan earthquake.
parameters used in the SG filter is the window size (half window for the filter) and a smooth polynomial degree (Chen et al., 2004, Diouf et al., 2015). Larger values of the window size produce a smoother result at the expense of flattening sharp peaks (Chen et al., 2004). A smaller value of the polynomial degree gives a smoother result but may introduce bias (Yang and Qi, 2017). However, if there are too many abnormal values due to frequent heavy cloud, the result of SG filter or improved SG filters (Hamunyela et al., 2013, 2016) could still have many valleys (Fig. 9). Taking the maximum values of every three consecutive bands (Max3 cloud filter) is another effective filter (Yang et al., 2013), but there still are abnormal values which cannot be eliminated sufficiently when cloud is frequent and heavy (Fig. 9). The Max3 cloud filter with outliers deletion proposed in this research can effectively eliminate small abnormal valleys and some stubborn
abnormal values and obtain an envelope curve to represent the original NDVI time series. To better diagnose the dynamic vegetation recovery after a major earthquake, appropriate cloud filters should be selected and tested in the study area to avoid potential bias or false judgement in the diagnosis results. 4.2. Effectiveness of DMI From the comparison of some specific locations using high-resolution satellite images, we found that both COI and Adjusted ED methods face some limitations in the diagnosis of vegetation recovery dynamics. The COI values can be easily affected by vegetation annual growth cycles. Although the COI method takes three consecutive bands as a time window to reduce the influence of phenological
57
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
changes, it may still generate false alarms when the vegetation growth is significantly delayed or advanced. The Adjusted ED method attempts to eliminate the influence of phenological changes by using a cross-correlagram algorithm, but it does not consider the accumulation of small differences induced by natural fluctuations. Therefore, the Adjusted ED method may exaggerate the interannual difference. The DMI proposed in this paper takes full advantages of the DTW and LBD algorithms to evaluate the vegetation recovery of each year after the earthquake. The use of DTW and LBD can effectively eliminate errors induced from phenological changes and the accumulation of small differences. The DMI created based on these two algorithms can better differentiate the annual differences in NDVI curves that are caused by real disturbances while avoiding false judgements on those differences due to normal phenological changes and natural fluctuations. It should be noted that the accuracy of the DMI relies on the data quality (or the selection of cloud filters) and the grid size (Wg) used in the calculation. The applicability of this method in other regions with different ecological characteristics requires further investigation. 4.3. Challenge in vegetation recovery after major earthquake The major earthquake occurred in mountainous areas usually cause wide spread landslides and severe disturbs to local vegetation cover and growth. In the following years after the earthquake, the locations with damaged and non-recovered vegetation may further
Fig. 6. Assessment of difference measurement trend.
Fig. 7. Vegetation change of Ma An Shi landslide area. 58
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 8. Vegetation change of typical forest area without disturbance.
series, and more effective and robust cloud filter methods are still needed. Different from many existing methods, such as those in references Lin et al. (2006, 2004, 2008, 2015), this method examines the vegetation state of a whole year period instead of some specific dates. The algorithm of DMI calculation considers a benchmark obtained from historical samples that were not influenced by disturbances (instead of specific year or date) and compares the investigated years after disturbances to the benchmark. The differences within the various investigated years are therefore comparable as they all calculated based on the same benchmark. The DTW-LBD algorithm proposed in this study overcomes some issues in existing methods. Using traditional Euclidean distance can generate misleading results because the interannual phenological change may amplify the difference in NDVI time series of different years. Using DTW-LBD algorithm allows annual NDVI time series to be stretched or compressed along the time axis to accommodate sequences that are similar but out of phase. This algorithm lowers the probability of false identification of difference between years. Moreover, we compared relatively regional vegetation state in various years with the historical benchmark and defined DMI. As a result, we diagnosed dynamic vegetation recovery and generated a series of interannual spatial distributions of regional vegetation conditions. The high temporal resolution of the MODIS products enables economical and efficient vegetation recovery assessment over a vast region for ecological monitoring after a major earthquake. Although this approach was developed to monitor vegetation recovery in mountainous regions affected by a major earthquake, it can be effectively used for vegetation analysis with other disturbances.
increase the likelihood of landslides occurrence (Cui et al., 2012), and this kind of impact due to post-seismic landslide activities may last for a long time, from years to decades (Yang et al., 2018). In the mountainous area, loose materials are very unstable after the major earthquake. With frequency and intensive rainfall events, new debris flow and landslides occurred and had negative influence on the nonlandslides areas. Moreover, land use types of many areas were changed after the earthquake. Some forest lands were changed to residential resettlement and farmlands, and new roads were constructed. Although the average vegetation condition at landslide sites has been recovering, the dispersion of vegetation condition is still large. The study area has been experiencing a fluctuating vegetation recovery condition, and more attention should be given to the vegetation monitoring and the prevention of subsequent hazards linked with vegetation degradation. 5. Conclusion The NDVI-based analysis often faces the issue of noises in signals particularly the effect of heavy and frequent cloud. The NDVI time series in the study area show frequent sharp drops that are difficult to be eliminated. The widely used SG filter or the Max3 filter may not be effective when the clouds are dense and often present. An improved practical filter is proposed to eliminate small abnormal valleys and stubborn abnormal values. Using Max3 cloud filter with outliers deletion, we can obtain the approximate envelopes of the original time series that can better reflect the actual vegetation status. It should be noted that the vegetation recovery assessment using DTW-LBD algorithm is highly sensitive to the process of the original NDVI time
59
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 9. Comparison of results using different filters.
Acknowledgments
Remote Sens. 2, 87–95. Chong, X., Fuchu, D., Jian, C., Xin-bin, T., Ling, X., Wei-chao, L., Wei, T., Yan-bo, C., Xin, Y., 2009. Identification and analysis of secondary geological hazards trig-gered by a magnitude 8.0 Wenchuan earthquake. J. Remote Sens. 13, 754–762. Cuevas-González, M., Gerard, F., Balzter, H., Riano, D., 2009. Analysing forest recovery after wildfire disturbance in boreal Siberia using remotely sensed vegetation indices. Glob. Change Biol. 15, 561–577. Cui, P., Lin, Y.-M., Chen, C., 2012. Destruction of vegetation due to geo-hazards and its environmental impacts in the Wenchuan earthquake areas. Ecol. Eng. 44, 61–69. Diouf, A.A., Brandt, M., Verger, A., Jarroudi, M.E., Djaby, B., Fensholt, R., Ndione, J.A., Tychon, B., 2015. Fodder biomass monitoring in sahelian rangelands using phenological metrics from FAPAR time series. Remote Sens. 7, 9122–9148. Ge, Y., Xu, J., Liu, Q., Yao, Y., Wang, R., 2009. Image interpretation and statistical analysis of vegetation damage caused by the Wenchuan earthquake and related secondary disasters. J. Appl. Remote Sens. 3 031660–031660 031619. Gouveia, C., DaCamara, C., Trigo, R., 2010. Post-fire vegetation recovery in Portugal based on spot/vegetation data. Nat. Hazards Earth Syst. Sci. 10, 673–684. Hamunyela, E., Verbesselt, J., Herold, M., 2016. Using spatial context to improve early detection of deforestation from Landsat time series. Remote Sens. Environ. 172, 126–138. Hamunyela, E., Verbesselt, J., Roerink, G., Herold, M., 2013. Trends in spring phenology
This research was supported by National Natural Science Foundation of China under grants 41671503 and 41505134. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ecolind.2018.06.026. References Chen, J., Jönsson, P., Tamura, M., Gu, Z., Matsushita, B., Eklundh, L., 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sens. Environ. 91, 332–344. Chen, Y.-R., Chen, J.-W., Hsieh, S.-C., Ni, P.-N., 2009. The application of remote sensing technology to the interpretation of land use for rainfall-induced landslides based on genetic algorithms and artificial neural networks. IEEE J. Sel. Top. Appl. Earth Obs.
60
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
259–273. Sun, M., 2012. Research on discords detect on time series based on distance and density. Comput. Eng. Appl 11–17+ 22. Tømmervik, H., Høgda, K.A., Solheim, I., 2003. Monitoring vegetation changes in Pasvik (Norway) and Pechenga in Kola Peninsula (Russia) using multitemporal Landsat MSS/TM data. Remote Sens. Environ. 85, 370–388. Turner, B.L., Lambin, E.F., Reenberg, A., 2007. The emergence of land change science for global environmental change and sustainability. Proc. Natl. Acad. Sci. 104, 20666–20671. Van Leeuwen, W.J., 2008. Monitoring the effects of forest restoration treatments on postfire vegetation recovery with MODIS multitemporal data. Sensors 8, 2017–2042. Wang, F., Zhou, L.-J., Liu, B., Zhang, W., Pan, F.-M., Luo, J.-G., 2009. Study on variation of forest landscape and indirect loss in Sichuan seriously disaster area by Wenchuan Earthquake. J. Soil Water Conserv. 5, 015. Wang, M., Liu, M., Yang, S., Shi, P., 2014a. Incorporating triggering and environmental factors in the analysis of earthquake-induced landslide hazards. Int. J. Disaster Risk Sci. 5, 125–135. Wang, M., Yang, W., Shi, P., Xu, C., Liu, L., 2014b. Diagnosis of vegetation recovery in mountainous regions after the Wenchuan earthquake. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 7, 3029–3037. Yang, W., Qi, W., 2017. Spatial-temporal dynamic monitoring of vegetation recovery after the Wenchuan earthquake. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 10, 868–876. Yang, W., Qi, W., Zhou, J., 2018. Decreased post-seismic landslides linked to vegetation recovery after the 2008 Wenchuan earthquake. Ecol. Ind. 89, 438–444. Yang, W., Wang, M., Shi, P., 2013. Using MODIS NDVI time series to identify geographic patterns of landslides in vegetated regions. IEEE Geosci. Remote Sens. Lett. 10, 707–710. Zhang, W., Lin, J., Peng, J., Lu, Q., 2010. Estimating Wenchuan Earthquake induced landslides based on remote sensing. Int. J. Remote Sens. 31, 3495–3508. Zhu, Y., Shasha, D., 2003. Warping indexes with envelope transforms for query by humming. Pages 181–192 in Proceedings of the 2003 ACM SIGMOD International Conference on Management of Data.
of western european deciduous? forests. Remote Sens. 5, 6159–6179. Jiang, W.G., Jia, K., Wu, J.J., Tang, Z.H., Wang, W.J., Liu, X.F., 2015. Evaluating the vegetation recovery in the damage area of Wenchuan earthquake using MODIS data. Remote Sens. 7, 8757–8778. Keogh, E., Ratanamahatana, C.A., 2005. Exact indexing of dynamic time warping. Knowl. Inf. Syst. 7, 358–386. Li, C., Wang, M., Liu, K., Xie, J., 2018. Topographic changes and their driving factors after 2008 Wenchuan earthquake. Geomorphology. Lin, C.Y., Lo, H.M., Chou, W.C., Lin, W.T., 2004. Vegetation recovery assessment at the Jou-Jou Mountain landslide area caused by the 921 Earthquake in Central Taiwan. Ecol. Model. 176, 75–81. Lin, W.-T., Lin, C.-Y., Chou, W.-C., 2006. Assessment of vegetation recovery and soil erosion at landslides caused by a catastrophic earthquake: a case study in Central Taiwan. Ecol. Eng. 28, 79–89. Lin, W.T., Chou, W.C., Lin, C.Y., 2008. Earthquake-induced landslide hazard and vegetation recovery assessment using remotely sensed data and a neural network-based classifier: a case study in central Taiwan. Nat. Haz. 47, 331–347. Liu, B., Fu-Zhong, W.U., Jian, Z., 2008. Key issues in restoration on earthquake-damaged ecosystem at the ecotone between dry valley and montane forest of the Minjiang River. Acta Ecol. Sin. Lu, T., Zeng, H., Luo, Y., Wang, Q., Shi, F., Sun, G., Wu, Y., Wu, N., 2012. Monitoring vegetation recovery after China’s May 2008 Wenchuan earthquake using Landsat TM time-series data: a case study in Mao County. Ecol. Res. 27, 955–966. Mantovani, F., Soeters, R., Van Westen, C., 1996. Remote sensing techniques for landslide studies and hazard zonation in Europe. Geomorphology 15, 213–225. Metternicht, G., Hurni, L., Gogu, R., 2005. Remote sensing of landslides: An analysis of the potential contribution to geo-spatial systems for hazard assessment in mountainous environments. Remote Sens. Environ. 98, 284–303. Ni, Z., He, Z., Zhao, Y., Gao, H., Cai, K.-K., Wang, L., 2009. Study on vegetation coverage changes in Dujiangyan before and after Wenchuan earthquake. Res. Soil Water Conserv. 16, 45–48. Röder, A., Hill, J., Duguy, B., Alloza, J.A., Vallejo, R., 2008. Using long time series of Landsat data to monitor fire events and post-fire dynamics and identify driving factors. A case study in the Ayora region (eastern Spain). Remote Sens. Environ. 112,
61
Contents lists available at ScienceDirect
Ecological Indicators journal homepage: www.elsevier.com/locate/ecolind
Using NDVI time series to diagnose vegetation recovery after major earthquake based on dynamic time warping and lower bound distance ⁎
Xuelei Zhanga, Ming Wanga,b, , Kai Liua,b, Jun Xieb, Hong Xua,
T
⁎
a State Key Laboratory of Earth Surface Processes and Resource Ecology/Academy of Disaster Reduction and Emergency Management, Faculty of Geographical Science, Beijing Normal University, Beijing 100875, China b Key Laboratory of Environmental Change and Natural Disasters, Ministry of Education, Beijing Normal University, Beijing 100875, China
A R T I C LE I N FO
A B S T R A C T
Keywords: Earthquake Vegetation recovery NDVI Dynamic time warping Lower bound distance Difference measurement index
Major earthquake occurred in mountainous areas usually cause large number of landslides that lead to severe impact to local vegetation cover and growth. The negative influence to vegetation may last for many year and vegetation recovery may experience dynamic fluctuation. Existing methods for vegetation recovery diagnosis face difficulty in capturing the dynamic behaviours both within and between years that makes the interannual comparison impossible. This paper proposes a new method to diagnose regional vegetation recovery after a major earthquake by defining a difference measurement index (DMI) using MODIS NDVI time series at 8-day interval. This differs from many existing methods in its quantification of the difference between the studied time series and historical samples, by using a proposed algorithm consisting of lower bound distance and dynamic time warping. This algorithm can better differentiate vegetation disturbance from its natural fluctuation. Second, the method investigates relatively regional vegetation recovery via a dynamic index, the DMI. Vegetation conditions in different years can be compared with a historical benchmark and measured by DMI. This makes it possible to diagnose dynamic vegetation recovery and generate a series of interannual spatial distributions of regional vegetation state.
1. Introduction When major earthquakes strike in mountainous regions, a large number of landslides, debris flows and quake lakes often occur and destroy vast areas of vegetation cover (Liu et al., 2008). After the Wenchuan Earthquake, several studies of regional ecological change and vegetation damage assessments were carried out using satellite data at multi-temporal scales (Ge et al., 2009). Damage to vegetation can be reflected in satellite images, for example via decreased Normalized Difference Vegetation Index (NDVI) values (Wang et al., 2014b). Much research has been undertaken concerning the immediate impact of earthquake-induced landslides on vegetation (Mantovani et al., 1996, Metternicht et al., 2005, Chen et al., 2009, Chong et al., 2009, Cui et al., 2012, Wang et al., 2014b). However, existing methods have difficulty in assessing the dynamics of vegetation recovery. More efforts should be paid to develop new methods for revealing the dynamic process of vegetation recovery and conduct case studies for crosscomparison of various methods. Remote sensing data can effectively diagnose vegetation recovery of mountainous regions (Tømmervik et al., 2003, Turner et al., 2007, Ni et al., 2009, Wang et al., 2009, Lu
et al., 2012), although researchers always face challenges from unsatisfactory data quality (e.g., cloud cover) and natural variability (e.g., changes in phenology) (Lu et al., 2012). The vegetation recovery rate (VRR) (Lin et al., 2006) method compares vegetation changes at a few cross-sectional data before and after a major disturbance. It provides an indicator of vegetation changes over time, however, it cannot give a complete picture of vegetation change since the occurrence of disturbance. The COI (Count of Outlier Index) (Wang et al., 2014b) determines an outlier if the NDVI value is double standard deviation less than past mean at a certain time window. The total number of outliers (the COI value) reflects the overall vegetation recovery status since the occurrence of the major disturbance, but it may generate false detection of outliers due to interannual phenological changes. Another approach based on the time series of satellite data using regression algorithms is to assess the trend of vegetation change and recovery rate (Röder et al., 2008, Van Leeuwen, 2008, Gouveia et al., 2010). However, the change in trends may significantly affected by the insufficient data length and the presence of abnormal values in original data. Very often, filters are used to eliminate the effect of cloud and abnormal points (e.g., cloud cover)
⁎ Corresponding authors at: State Key Laboratory of Earth Surface Processes and Resource Ecology/Academy of Disaster Reduction and Emergency Management, Faculty of Geographical Science, Beijing Normal University, Beijing 100875, China (M. Wang for general correspondence and H. Xu for data). E-mail addresses: [email protected] (M. Wang), [email protected] (H. Xu).
https://doi.org/10.1016/j.ecolind.2018.06.026 Received 23 October 2017; Received in revised form 6 June 2018; Accepted 12 June 2018 1470-160X/ © 2018 Elsevier Ltd. All rights reserved.
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 1. Geographic location of study area (eastern Pingwu County).
(Cuevas-González et al., 2009). Savitzky–Golay filter (SG filter) is widely used based on a weighted moving average algorithm (Chen et al., 2004), but the determination of filter parameters usually influences analysis results (Chen et al., 2004, Diouf et al., 2015, Yang and Qi, 2017). Taking the maximum values of every three consecutive bands (Max3 cloud filter) is another effective filter (Yang et al., 2013), but frequent and heavy cloud may not be eliminated effectively. In processing time series data, we need to eliminate obvious abnormal signals (e.g., cloud cover) while retaining signals of sudden change that may be induced by vegetation damage. Therefore, careful processing is required. In this paper, we develop a new method based on 8-day interval MODIS (Moderate Resolution Imaging Spectroradiometer) NDVI time series to diagnose vegetation recovery after a major earthquake. This differs from many existing methods in the way of quantifying the difference between a time series and historical samples by using a proposed DTW-LBD algorithm that consists of LBD (lower bound distance) and DTW (dynamic time warping). LBD is used for increasing the substantial difference between time series while neglecting small distinctions (i.e., NDVI variation caused by moderately increased or reduced seasonal rainfall), and time warping (TW) similarity measure that allows sequences to be stretched or compressed along the time axis. The DTW algorithm (Keogh and Ratanamahatana, 2005) allows elastic shifting of the time axis to accommodate sequences that are similar but out of phase. DTW is used for eliminating any errors induced by interannual phenological change (i.e., growing season delay due to extended winter period). The DTW-LBD algorithm therefore can better differentiate vegetation disturbance from its natural fluctuation and enables diagnosis of dynamic vegetation recovery and generation of a series of interannual spatial distributions of regional vegetation state. The objectives of this paper are to, first, develop a new method that can assess the dynamic change of vegetation recovery after major earthquake by minimizing the influence of interannual phenological changes and persistent small difference, and second, use a case study of the Wenchuan earthquake stricken area to verify the validity and applicability of the proposed method.
2. Materials and methods 2.1. Study area The Wenchuan region in southwest China is ecologically diverse and but fragile. At 2:28p.m. local time on May 12th, 2008, a magnitude 8.0 earthquake occurred in Sichuan Province, followed by thousands of associated aftershocks, resulting in more than 200,000 landslides over an area of 110,000 km2 (Wang et al., 2014b). Eastern Pingwu County was one of the most severely affected areas during the Wenchuan Earthquake, with Modified Mercalli Intensity scale IX and X. Southeast monsoon dominants the local climate. Most rainfall concentrate from July to October with annual average precipitation of 800 mm and a standard deviation of 167 mm (Wang et al., 2014a, Li et al., 2018). The average temperature is 14.8 °C with daily minimum of −4 °C and maximum of 39.1 °C. This geographic extent shown in Fig. 1 was selected to test the proposed DTW-LBD algorithm for vegetation recovery diagnosis. 2.2. Materials The MODIS sensor is on board the Terra and Aqua satellites. We used MODIS MOD09Q1 products (Wang et al., 2014b), with spatial resolution 250 meters and temporal resolution nearly 8 days. MOD09Q1 contains both red light (red) and near-infrared (NIR) bands. We converted the projection to WGS-84, UTM48N according to the location of the study area (Wang et al., 2014b). In that area, 690 images were acquired from January 1, 2001 to December 31, 2015. The NDVIs for all cells were calculated based on reflectance in the NIR and red bands (Zhang et al., 2010, Wang et al., 2014b)
NDVI = (NIR−red )/(NIR + red )
(1)
Cloud, aerosol, and vapor in the atmosphere usually cause an abnormal decrease of the NDVI, thereby not reflecting the actual vegetation condition. Residues that were not detected during the masking process of atmospheric noise also generate abnormal valleys in the time
53
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
series. Here we use Max3 cloud filter with outliers deletion to achieve this goal by taking the following steps. Firstly, per reference (Yang et al., 2013), we took the maximum value of every three consecutive bands to form a new time series (Max3 cloud filter):
NDVIi = max (NDVIi − 1, NDVIi , NDVIi + 1, )
(2)
In this formula, i means the ith time band in original NDVI time series. Secondly, we remove the NDVI values higher than 0.99, and use the average of the left and right NDVI values to substitute the abnormal value. Thirdly, we define an abnormal valley when the NDVI experiences a sudden reduction followed by an immediately sudden recovery. If the value of the previous NDVI is 0.2 higher than the current point while the value of the next NDVI is also 0.2 higher, we define this point as an abnormal valley. In addition, there is a phenomenon that abnormal valley may have more than one point. We also use the average of the left and right NDVI of the abnormal valley to substitute the values in the valley.
Fig. 2. Typical cases of LBD calculation.
minimum difference between two sequences, DTW maps each element of a sequence to one or more neighbouring elements of another sequence. Given any two non-null sequences A and B , the time warping LBD is defined as follows.
2.3. Grid representation and LBD The time series data can be divided by predetermined grid cells and transformed into a form of grid representation (Sun, 2012). A = {a1, a2 , ⋯, an} is the original time series and ai is the value at time i . If the grid width (Wg) is given, the grid representation of A can be defined as GRA = {〈AR1 , ARB1 〉, 〈AR2 , ARB2 〉, ⋯, 〈ARn , ARBn 〉} , in which ARi is the row of ith value of A in the grid and ARBi represents which part of a grid into which the ith value is entered. ABi is an integer denoting the rower. If the ith value of A is above the center line of a grid, ARBi = 1; otherwise ARBi = 0 . The centre line of a grid is calculated by the average of upper and lower bound. In contrast with the widely used Euclidean distance, LBD is used to quantify the difference between two time series A and B based on their grid representations GRA = {〈AR1 , ARB1 〉, 〈AR2 , ARB2 〉, ⋯, 〈ARn , ARBn 〉} and GRB = {〈BR1, BRB1 〉, 〈BR2 , BRB2 〉, ⋯, 〈BRn , BRBn 〉} .
DTWLBD (A, B ) =
Here, First (A) denotes the ‘current’ first element in A, and Rest (A) denotes the sequence consisting of the elements in A other than First (A) . This equation shows an iteration procedure in which the ‘min’ term repeats the DTW–LBD algorithm and the First ( ) and Rest ( ) functions keep updating the first and remaining sequences (Zhu and Shasha, 2003).
n
2.5. Measurement of vegetation recovery
n
⎨ ∑ (AR −BR −1)2Wg 2 , ifAR ≠ BR , i i i i ⎪ ⎪ i=1 ⎪ and (ARi −BRi )(ARBi −BRBi ) ⩽ 0 ⎪ ⎩ 0, if ARi = BRi
2 ⎧ DTWLBD (A, Rest (B )) First (B )) + min DTWLBD 2 (Rest (A), B ) ⎨ 2 ⎩ DTWLBD (Rest (A), Rest (B ))
(4)
⎧ 2 2 ⎪ ∑ (ARi −BRi ) Wg , ifARi ≠ BRi , ⎪ i=1 ⎪ ⎪ and (ARi −BRi )(ARBi −BRBi ) > 0 LBD (A, B ) =
LBD 2 (First (A),
Based on the proposed DTW–LBD algorithm, a difference measurement index (DMI) can be defined: ym + x
DMIx (k ) =
(3)
Fig. 2 uses five typical cases to illustrate the rules in formula (3). When A and B is on the same grid as showed (case (a) and (b)), the LBD of A and B is zero. When A and B is on different grids, there are three typical cases in which the difference between A and B is approximately estimated by counting the grids between A and B. Using the LBD algorithm, the difference of the ith observed values is decreased. In case (c) and (d), the LBD is estimated to be one grid (Wg) as (ARBi −BRBi ) ⩽ 0 , while the LBD is two grids (2 Wg) for case (e) as (ARBi −BRBi ) > 0 .
∑y = y1 DTWLBD (Tyi, Tyj ) Cm2 + 1
ym
−
∑y = y1 DTWLBD (Tyi, Tyj ) Cm2
(5)
Here, DMIx (k ) denotes the vegetation recovery status at location k ym in year x . The term ∑y = y1 DTWLBD (Tyi, Tyj )/Cm2 represents the average DTW–LBD of any two time series Tyi and Tyj in samples from year y1 to ym . This is considered a benchmark of DTW–LBD obtained from historical samples that were not influenced by disturbances (e.g., earthquakes). In this article, y1 = 2001 and ym = 2007 . The term ym + x ∑y = y1 DTWLBD (Tyi, Tyj )/Cm2 + 1 denotes the average DTW–LBD of any two time series Tyi and Tyj in samples from year y1 to ym plus year x . In this article, x = 2008, 2009, ⋯, or 2015. This term differs from the benchmark by adding one more time series of the investigated year. Therefore, the difference between these two terms is larger if the added time series from year x is significantly different from the historical average. DMIx (k ) takes the DTW-LBD difference comparing to historical average. This makes DMIx (k ) comparable between years. As a result, the bigger DMI is, the worse vegetation recovered. Linear regression can be done to analyze the difference measurement trend (DMT) based on calculated annual DMIs, as follows.
2.4. DTW of time series with LBD There is always natural interannual fluctuation of vegetation. Time series patterns may change and lead to dramatic differences if calculated using simply the distance at the same time. Here, we use a time warping (TW) similarity measure that allows sequences to be stretched or compressed along the time axis. The DTW algorithm (Keogh and Ratanamahatana, 2005) allows elastic shifting of the time axis to accommodate sequences that are similar but out of phase. To find the
DMIx (k ) = DMT × x + β0 54
(6)
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Step 1: Prepare MODIS images and use the NIR and red bands to calculate the NDVI time series. Step 2: Apply the filtering algorithm to eliminate possible cloud effects. Step 3: Determine a grid size and transform the NDVI time series into grid format, including rows and 0/1, denoting locations on a grid. Step 4: Calculate DMI for all locations based on the historical benchmark samples and samples from the investigated year, using the DTW–LBD algorithms. Step 5: Calculate DMT by using the linear regression based on DMI results. Step 6: Perform DMI and DMT mapping to illustrate vegetation status and recovery trends.
Table 1 High-resolution satellite image in study area. Sensor
Date
Resolution (panchromatic)
IKONOS
November 2, 2002 and November 8, 2002 April 23, 2010 and August 15, 2010 January 29, 2014 August 20, 2015
1m
WV02+qb ZY-3 GF-2
0.5 m 2.5 m 0.8 m
Here, β0 and β1 are the regression coefficients. If DMT < 0 , the DMI at the analysed location decreases with the increase of years and vegetation tends to be more similar to the historical benchmark, implying a positive trend of recovery. If DMT > 0 , the vegetation tends to have a more significant difference relative to the historical benchmark, implying a negative trend of recovery. Therefore, the distribution of regional recovery trends can be obtained by using values of DMT , and it can be used for qualitative assessment of the recovery trend.
3. Results 3.1. Effect of DTW-LBD We chose some typical NVDI time series from the study area to demonstrate the effectiveness of the proposed method and compare with other methods. As shown in Fig. 3, the red line is one of the NDVI time series of 2001, and the blue line is one of the NDVI time series of 2014. The original NDVI time series have big difference due to the influence of cloud, and after application of the cloud filter, the difference is mainly related to the interannual phenological change. In term of vegetation condition, these two NDVI curves in 2001 and 2014 are considered to be similar as no significant difference is observed other than the phenological change. If we use difference methods to calculate the distance between these two NDVI time series, the results can be very different. As shown in Table 2, Euclid distance is highest. The adjusted Euclid distance (Yang and Qi, 2017) (Adjusted ED) applied the SG filter and a cross-correlogram spectral matching algorithm to eliminate the influences of cloud and interannual phenological changes, but still it is much higher than the
2.6. Accuracy assessment For assessing accuracy of this new method, we used some high-resolution satellite images (Table 1) and statistics analysis to verify the assessment results. We also plotted the NVDI distributions and did statistics analysis of their mean and variance of those picked areas in different years. Z-test and F-test are used to check whether significant difference exists in the mean and variance of NDVI among different years, respectively. 2.7. Step-by-step procedure The following steps were used to generate maps of vegetation recovery in various years, based on the MODIS NDVI time series.
Fig. 3. Influence of interannual phenological change.
Table 2 Comparison of distance calculated using different methods.
Case with phenological change (Fig. 3) Case with persistent small difference (Fig. 4)
Euclid Dist.
Adjusted ED
Dist.w/DTW
DTWLBD (Wg = 0.01)
DTWLBD (Wg = 0.05)
DTWLBD (Wg = 0.1)
1.8789 1.5803
1.1012 0.6003
0.4547 0.1784
0.4144 0.1371
0.3122 0.0500
0.2236 0
55
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 4. Influence of persistent small difference.
albeit with some fluctuation. Many non-landslide areas showed a negative recovery trend, indicating that earthquake-stricken areas still faced major challenges because of subsequent geo-hazards or man-made disturbances.
proposed distance with DTW and LBD. If we use NDVI time series in 2001 as the base pre-seismic NDVI time series and the distance between these two NDVI time series to represent vegetation recovery, the Euclid distance and Adjusted ED may not provide an accurate assessment. So DTW could eliminate the influence of interannual phonological changes. Persistent small difference between two time series can also be accumulated into a big difference. As shown in Fig. 4, the red line is one NDVI time series of 2003, and the blue line is one NDVI time series of 2006. There is persistent small difference between these two NDVI time series, and obviously the vegetation condition in these two years is considered to be similar. As shown in Table 2, Euclid distance is highest. With the use of LBD algorithm, this kind of persistent small difference can be eliminated.
3.3. Accuracy assessment The DTW-LBD method is compared with the COI method proposed by Wang et al., (2014b) and the adjusted ED method proposed by Yang and Qi (2017). We picked two areas for the purpose of verification: Ma An Shi landslide (I in Fig. 6) and a typical forest without disturbance (VII in Fig. 6). From the high-resolution satellite images of 2010 in Ma An Shi landslide (Fig. 7), we can see that the landslides caused vegetation deterioration. The vegetation was recovering gradually in the following years such as 2014 and 2015. The mean of NDVI is increasing. In 2015, the mean NDVI is even higher than that in 2002 (significant at p < 0.05), but the variance is also higher (significance at p < 0.05), implying that many local areas experienced vegetation degradation although the overall vegetation was recovering. The results of DMI well demonstrated this dynamic recovery process. The DMI results were compared with the COI and adjusted ED results, suggesting that the DMI results in general agreed with the COI results while the COI results may sometimes exaggerate the difference. The adjusted ED results did not provide sufficient differentiation in terms of spatial patterns. We also selected an area of typical forest that experienced no obvious disturbance (Fig. 8). The mean and variance of the NDVI in 2014 and 2015 remain similar (significance at p < 0.05), which agreed with the DMI results. However, the other two methods exaggerated the difference in some years, which may not accurately reflect the actual situations.
3.2. Vegetation recovery diagnosis We used 690 images in the study area to create the NDVI time series, spanning 2001–2015. By taking samples from 2001 to 2007 as the benchmark, annual DMIs were obtained for 2008–2015. Fig. 5 shows that the study area had dynamic change of vegetation recovery after the major earthquake in 2008. DMIs in all eight figures are at the same scale, so interannual differences are comparable. Large DMI values in 2008 correspond to locations of large landslides, mainly triggered by the earthquake. That means severe vegetation disturbances right after the major earthquake concentrated in the areas with large co-seismic landslides. The heaviest damage to vegetation (in red or orange colour) occurred in the northern part of the area, where two top ten largest landslides are located. These include Ma An Shi (west) and Zheng Jia Shan (east). Starting from 2009, the vegetation in those landslides areas recovered gradually, while vegetation disturbance started to appear in those non-landslide areas, particularly in the southern part of the study area. The vegetation recovery in the study area shows an overall fluctuating and spatial shift. Fig. 6 shows the recovery trend of the study area for the seven years after the Wenchuan earthquake. All large-scale landslides in the study area (I-VI circled in green) had a positive recovery trend, and vegetation of Shi Ban Po landslide, Wen Jia Ba landslide, Zheng Jia Wan landslide recovered (significance at p < 0.1). These landslide areas had heavy vegetation damage in 2008 but gradually recovered their vegetation,
4. Discussion 4.1. Effect of cloud filter selection Savitzky–Golay filter (SG filter) is a widely used method to resolve cloud noise in NDVI time series using a weighted moving average filter with the weighting given as a polynomial (Chen et al., 2004). The key
56
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 5. Assessment of annual vegetation change after Wenchuan earthquake.
parameters used in the SG filter is the window size (half window for the filter) and a smooth polynomial degree (Chen et al., 2004, Diouf et al., 2015). Larger values of the window size produce a smoother result at the expense of flattening sharp peaks (Chen et al., 2004). A smaller value of the polynomial degree gives a smoother result but may introduce bias (Yang and Qi, 2017). However, if there are too many abnormal values due to frequent heavy cloud, the result of SG filter or improved SG filters (Hamunyela et al., 2013, 2016) could still have many valleys (Fig. 9). Taking the maximum values of every three consecutive bands (Max3 cloud filter) is another effective filter (Yang et al., 2013), but there still are abnormal values which cannot be eliminated sufficiently when cloud is frequent and heavy (Fig. 9). The Max3 cloud filter with outliers deletion proposed in this research can effectively eliminate small abnormal valleys and some stubborn
abnormal values and obtain an envelope curve to represent the original NDVI time series. To better diagnose the dynamic vegetation recovery after a major earthquake, appropriate cloud filters should be selected and tested in the study area to avoid potential bias or false judgement in the diagnosis results. 4.2. Effectiveness of DMI From the comparison of some specific locations using high-resolution satellite images, we found that both COI and Adjusted ED methods face some limitations in the diagnosis of vegetation recovery dynamics. The COI values can be easily affected by vegetation annual growth cycles. Although the COI method takes three consecutive bands as a time window to reduce the influence of phenological
57
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
changes, it may still generate false alarms when the vegetation growth is significantly delayed or advanced. The Adjusted ED method attempts to eliminate the influence of phenological changes by using a cross-correlagram algorithm, but it does not consider the accumulation of small differences induced by natural fluctuations. Therefore, the Adjusted ED method may exaggerate the interannual difference. The DMI proposed in this paper takes full advantages of the DTW and LBD algorithms to evaluate the vegetation recovery of each year after the earthquake. The use of DTW and LBD can effectively eliminate errors induced from phenological changes and the accumulation of small differences. The DMI created based on these two algorithms can better differentiate the annual differences in NDVI curves that are caused by real disturbances while avoiding false judgements on those differences due to normal phenological changes and natural fluctuations. It should be noted that the accuracy of the DMI relies on the data quality (or the selection of cloud filters) and the grid size (Wg) used in the calculation. The applicability of this method in other regions with different ecological characteristics requires further investigation. 4.3. Challenge in vegetation recovery after major earthquake The major earthquake occurred in mountainous areas usually cause wide spread landslides and severe disturbs to local vegetation cover and growth. In the following years after the earthquake, the locations with damaged and non-recovered vegetation may further
Fig. 6. Assessment of difference measurement trend.
Fig. 7. Vegetation change of Ma An Shi landslide area. 58
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 8. Vegetation change of typical forest area without disturbance.
series, and more effective and robust cloud filter methods are still needed. Different from many existing methods, such as those in references Lin et al. (2006, 2004, 2008, 2015), this method examines the vegetation state of a whole year period instead of some specific dates. The algorithm of DMI calculation considers a benchmark obtained from historical samples that were not influenced by disturbances (instead of specific year or date) and compares the investigated years after disturbances to the benchmark. The differences within the various investigated years are therefore comparable as they all calculated based on the same benchmark. The DTW-LBD algorithm proposed in this study overcomes some issues in existing methods. Using traditional Euclidean distance can generate misleading results because the interannual phenological change may amplify the difference in NDVI time series of different years. Using DTW-LBD algorithm allows annual NDVI time series to be stretched or compressed along the time axis to accommodate sequences that are similar but out of phase. This algorithm lowers the probability of false identification of difference between years. Moreover, we compared relatively regional vegetation state in various years with the historical benchmark and defined DMI. As a result, we diagnosed dynamic vegetation recovery and generated a series of interannual spatial distributions of regional vegetation conditions. The high temporal resolution of the MODIS products enables economical and efficient vegetation recovery assessment over a vast region for ecological monitoring after a major earthquake. Although this approach was developed to monitor vegetation recovery in mountainous regions affected by a major earthquake, it can be effectively used for vegetation analysis with other disturbances.
increase the likelihood of landslides occurrence (Cui et al., 2012), and this kind of impact due to post-seismic landslide activities may last for a long time, from years to decades (Yang et al., 2018). In the mountainous area, loose materials are very unstable after the major earthquake. With frequency and intensive rainfall events, new debris flow and landslides occurred and had negative influence on the nonlandslides areas. Moreover, land use types of many areas were changed after the earthquake. Some forest lands were changed to residential resettlement and farmlands, and new roads were constructed. Although the average vegetation condition at landslide sites has been recovering, the dispersion of vegetation condition is still large. The study area has been experiencing a fluctuating vegetation recovery condition, and more attention should be given to the vegetation monitoring and the prevention of subsequent hazards linked with vegetation degradation. 5. Conclusion The NDVI-based analysis often faces the issue of noises in signals particularly the effect of heavy and frequent cloud. The NDVI time series in the study area show frequent sharp drops that are difficult to be eliminated. The widely used SG filter or the Max3 filter may not be effective when the clouds are dense and often present. An improved practical filter is proposed to eliminate small abnormal valleys and stubborn abnormal values. Using Max3 cloud filter with outliers deletion, we can obtain the approximate envelopes of the original time series that can better reflect the actual vegetation status. It should be noted that the vegetation recovery assessment using DTW-LBD algorithm is highly sensitive to the process of the original NDVI time
59
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
Fig. 9. Comparison of results using different filters.
Acknowledgments
Remote Sens. 2, 87–95. Chong, X., Fuchu, D., Jian, C., Xin-bin, T., Ling, X., Wei-chao, L., Wei, T., Yan-bo, C., Xin, Y., 2009. Identification and analysis of secondary geological hazards trig-gered by a magnitude 8.0 Wenchuan earthquake. J. Remote Sens. 13, 754–762. Cuevas-González, M., Gerard, F., Balzter, H., Riano, D., 2009. Analysing forest recovery after wildfire disturbance in boreal Siberia using remotely sensed vegetation indices. Glob. Change Biol. 15, 561–577. Cui, P., Lin, Y.-M., Chen, C., 2012. Destruction of vegetation due to geo-hazards and its environmental impacts in the Wenchuan earthquake areas. Ecol. Eng. 44, 61–69. Diouf, A.A., Brandt, M., Verger, A., Jarroudi, M.E., Djaby, B., Fensholt, R., Ndione, J.A., Tychon, B., 2015. Fodder biomass monitoring in sahelian rangelands using phenological metrics from FAPAR time series. Remote Sens. 7, 9122–9148. Ge, Y., Xu, J., Liu, Q., Yao, Y., Wang, R., 2009. Image interpretation and statistical analysis of vegetation damage caused by the Wenchuan earthquake and related secondary disasters. J. Appl. Remote Sens. 3 031660–031660 031619. Gouveia, C., DaCamara, C., Trigo, R., 2010. Post-fire vegetation recovery in Portugal based on spot/vegetation data. Nat. Hazards Earth Syst. Sci. 10, 673–684. Hamunyela, E., Verbesselt, J., Herold, M., 2016. Using spatial context to improve early detection of deforestation from Landsat time series. Remote Sens. Environ. 172, 126–138. Hamunyela, E., Verbesselt, J., Roerink, G., Herold, M., 2013. Trends in spring phenology
This research was supported by National Natural Science Foundation of China under grants 41671503 and 41505134. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ecolind.2018.06.026. References Chen, J., Jönsson, P., Tamura, M., Gu, Z., Matsushita, B., Eklundh, L., 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sens. Environ. 91, 332–344. Chen, Y.-R., Chen, J.-W., Hsieh, S.-C., Ni, P.-N., 2009. The application of remote sensing technology to the interpretation of land use for rainfall-induced landslides based on genetic algorithms and artificial neural networks. IEEE J. Sel. Top. Appl. Earth Obs.
60
Ecological Indicators 94 (2018) 52–61
X. Zhang et al.
259–273. Sun, M., 2012. Research on discords detect on time series based on distance and density. Comput. Eng. Appl 11–17+ 22. Tømmervik, H., Høgda, K.A., Solheim, I., 2003. Monitoring vegetation changes in Pasvik (Norway) and Pechenga in Kola Peninsula (Russia) using multitemporal Landsat MSS/TM data. Remote Sens. Environ. 85, 370–388. Turner, B.L., Lambin, E.F., Reenberg, A., 2007. The emergence of land change science for global environmental change and sustainability. Proc. Natl. Acad. Sci. 104, 20666–20671. Van Leeuwen, W.J., 2008. Monitoring the effects of forest restoration treatments on postfire vegetation recovery with MODIS multitemporal data. Sensors 8, 2017–2042. Wang, F., Zhou, L.-J., Liu, B., Zhang, W., Pan, F.-M., Luo, J.-G., 2009. Study on variation of forest landscape and indirect loss in Sichuan seriously disaster area by Wenchuan Earthquake. J. Soil Water Conserv. 5, 015. Wang, M., Liu, M., Yang, S., Shi, P., 2014a. Incorporating triggering and environmental factors in the analysis of earthquake-induced landslide hazards. Int. J. Disaster Risk Sci. 5, 125–135. Wang, M., Yang, W., Shi, P., Xu, C., Liu, L., 2014b. Diagnosis of vegetation recovery in mountainous regions after the Wenchuan earthquake. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 7, 3029–3037. Yang, W., Qi, W., 2017. Spatial-temporal dynamic monitoring of vegetation recovery after the Wenchuan earthquake. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 10, 868–876. Yang, W., Qi, W., Zhou, J., 2018. Decreased post-seismic landslides linked to vegetation recovery after the 2008 Wenchuan earthquake. Ecol. Ind. 89, 438–444. Yang, W., Wang, M., Shi, P., 2013. Using MODIS NDVI time series to identify geographic patterns of landslides in vegetated regions. IEEE Geosci. Remote Sens. Lett. 10, 707–710. Zhang, W., Lin, J., Peng, J., Lu, Q., 2010. Estimating Wenchuan Earthquake induced landslides based on remote sensing. Int. J. Remote Sens. 31, 3495–3508. Zhu, Y., Shasha, D., 2003. Warping indexes with envelope transforms for query by humming. Pages 181–192 in Proceedings of the 2003 ACM SIGMOD International Conference on Management of Data.
of western european deciduous? forests. Remote Sens. 5, 6159–6179. Jiang, W.G., Jia, K., Wu, J.J., Tang, Z.H., Wang, W.J., Liu, X.F., 2015. Evaluating the vegetation recovery in the damage area of Wenchuan earthquake using MODIS data. Remote Sens. 7, 8757–8778. Keogh, E., Ratanamahatana, C.A., 2005. Exact indexing of dynamic time warping. Knowl. Inf. Syst. 7, 358–386. Li, C., Wang, M., Liu, K., Xie, J., 2018. Topographic changes and their driving factors after 2008 Wenchuan earthquake. Geomorphology. Lin, C.Y., Lo, H.M., Chou, W.C., Lin, W.T., 2004. Vegetation recovery assessment at the Jou-Jou Mountain landslide area caused by the 921 Earthquake in Central Taiwan. Ecol. Model. 176, 75–81. Lin, W.-T., Lin, C.-Y., Chou, W.-C., 2006. Assessment of vegetation recovery and soil erosion at landslides caused by a catastrophic earthquake: a case study in Central Taiwan. Ecol. Eng. 28, 79–89. Lin, W.T., Chou, W.C., Lin, C.Y., 2008. Earthquake-induced landslide hazard and vegetation recovery assessment using remotely sensed data and a neural network-based classifier: a case study in central Taiwan. Nat. Haz. 47, 331–347. Liu, B., Fu-Zhong, W.U., Jian, Z., 2008. Key issues in restoration on earthquake-damaged ecosystem at the ecotone between dry valley and montane forest of the Minjiang River. Acta Ecol. Sin. Lu, T., Zeng, H., Luo, Y., Wang, Q., Shi, F., Sun, G., Wu, Y., Wu, N., 2012. Monitoring vegetation recovery after China’s May 2008 Wenchuan earthquake using Landsat TM time-series data: a case study in Mao County. Ecol. Res. 27, 955–966. Mantovani, F., Soeters, R., Van Westen, C., 1996. Remote sensing techniques for landslide studies and hazard zonation in Europe. Geomorphology 15, 213–225. Metternicht, G., Hurni, L., Gogu, R., 2005. Remote sensing of landslides: An analysis of the potential contribution to geo-spatial systems for hazard assessment in mountainous environments. Remote Sens. Environ. 98, 284–303. Ni, Z., He, Z., Zhao, Y., Gao, H., Cai, K.-K., Wang, L., 2009. Study on vegetation coverage changes in Dujiangyan before and after Wenchuan earthquake. Res. Soil Water Conserv. 16, 45–48. Röder, A., Hill, J., Duguy, B., Alloza, J.A., Vallejo, R., 2008. Using long time series of Landsat data to monitor fire events and post-fire dynamics and identify driving factors. A case study in the Ayora region (eastern Spain). Remote Sens. Environ. 112,
61