Prediction of Vegetation Change by Discrete Wavelet Decomposition Based on Remote Sensing Time Series Images

Prediction of Vegetation Change by Discrete Wavelet Decomposition Based on Remote Sensing Time Series Images

Yuehong Long Jianxin Qin* Ke Wang Yun Xue Ling Wang

School of Geographical Sciences, Hunan Normal University, Changsha 410000, China

College of Municipal and Geomatics Engineering, Hunan City University, Yiyang 413000, China

College of Information and Electronic Engineering, Hunan City University, Yiyang 413000, China

College of Economic and Trade, University of Electronic Science and Technology of China, Zhongshan Institute, Zhongshan 528499, China

Corresponding Author Email: 
12495@hunnu.edu.cn
Page: 
123-132
|
DOI: 
https://doi.org/10.18280/ts.400111
Received: 
9 November 2022
|
Revised: 
16 January 2023
|
Accepted: 
26 January 2023
|
Available online: 
28 February 2023
| Citation

© 2023 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).

OPEN ACCESS

Abstract: 

The development of remote sensing technology has accumulated a large number of remote sensing image time series data for human monitoring of surface vegetation change, which provides a basis for vegetation change prediction. In order to improve the prediction accuracy of vegetation change, this paper uses discrete wavelet to decompose remote sensing image sequences at multiple scales, to explore the difference of influence of different temporal scale change characteristics on vegetation spatio-temporal change prediction, and find the best decomposition scale for vegetation change prediction. In this paper, the research object is the MODIS 13Q1 EVI image data of Hunan Province from 2001 to 2021. The discrete wavelet is adopted to obtain multi-scale vegetation trend components and detailed component sequences, and then complete the LSTM modeling prediction and comparison. The following are the experimental findings: the predictive ability of the discrete wavelet decomposition sequence group is better than that of the original EVI time series to varying degrees. The order of prediction accuracy is: monthly scale > seasonal scale > annual scale > original EVI time series. Thus, it is of reference significance to the research of application scenarios of change prediction of other regionalized variables with multi-scale characteristics.

Keywords: 

discrete wavelet, remote sensing time Series images, multi-scale, spatio-temporal change prediction, LSTM

1. Introduction

The development of remote sensing technology has accumulated a large number of remote sensing image time series data for human monitoring of surface vegetation change, providing a reliable data source for vegetation change prediction [1, 2]. Affected by natural factors and human activities, vegetation change exhibits long-term, seasonal periodicity and short-term multi-scale changes characteristics [3-6]. The long-term remote sensing image data obtained by remote sensing technology records the complex change process of the vegetation on the earth’s surface. However, the variation features of different scales contained in the image sequences make the data features fuzzy due to mutual interference, which leads to the low accuracy of vegetation prediction to a certain extent. Therefore, the separation of these different scale features and the exploration of their applicability to vegetation change prediction are of great significance to the study of vegetation change prediction.

With the development of wavelet theory in recent years, wavelet transform (WT) has been widely used in signal processing, image segmentation, and time-series data analysis due to its multiresolution analysis capability [7-9]. The discrete wavelet has time series data multiresolution decomposition and reconstruction capabilities. Time series data can be decomposed into time series data groups representing different time scales through wavelet transform. At present, some studies have shown that the prediction of time series data transformed by discrete wavelet can obtain relatively high accuracy [10-12]. These decomposed sequence group data can better retain the change trend features and detail features on different time scales, and enhance the prediction model’s ability to understand high-pass detail information and low-pass trend information at different scales. According to the domestic and foreign research on vegetation change, the role of wavelet transform’s multi-scale analysis theory in monitoring the dynamic changes of vegetation is demonstrated. It can be found that the wavelet analysis of a set of vegetation index time series can mine the dynamic change characteristics of vegetation within a year, as well as the dynamic change trend of vegetation between years [3, 13-15]. With the decomposition and reconstruction ability of vegetation time series data, discrete wavelet with great potential in vegetation prediction can provide a basis for multi-scale predicting of vegetation changes and analyzing vegetation change prediction capabilities of different scales of vegetation change features.

Early research on vegetation change prediction mainly established a series of simulation models by studying the relationship between global climate change and vegetation, and then related research used climate and other related data to simulate and predict vegetation. With the development of the spatio-temporal series prediction method, especially after the introduction of the deep learning mechanism, the deep learning method and its data-driven approach have shown relatively strong nonlinear learning capabilities. In this case, the spatio-temporal series prediction method can better improve the autocorrelation of explanatory variables and reduce uncertainty, which can achieve good prediction results in many fields such as finance, meteorology, and transportation. In particular, the long short term memory networks model (LSTM) has excellent performance in time series prediction [16], and has become a research hotspot in the field of deep learning in recent years. The model maintains a good memory for time series with a long-time span, and can solve the gradient disappearance and gradient explosion problems existing in recurrent neural networks [17, 18]. At present, the model has active performance in the prediction research of social and economic fields such as stock price index [19], exchange rate [20], transportation [21], energy consumption [22], and agricultural product price [23]. In the field of natural sciences, LSTM is used in the prediction research of water resources [24, 25], lightning [26], and air pollution [27], which has achieved excellent results. Some scholars have introduced LSTM into the research of vegetation change prediction [28, 29]. However, LSTM is still rare in the study of vegetation change, especially in the study of temporal multi-scale decomposition sequence vegetation change prediction.

This study uses the MODIS 13Q1 EVI time series images of Hunan Province from 2001 to 2021. The role of decomposition sequences in vegetation change prediction is verified by combining prediction of multi-scale decomposition sequences of vegetation time series. The differences in the applicability of different time-scale change features to vegetation change prediction are studied. The decomposition sequence group with the best prediction performance will be found, so that the vegetation change prediction results can achieve higher accuracy.

2. Data and Preprocessing

2.1 Vegetation remote sensing data

In this paper, the vegetation change prediction research uses the 2001-2021 MOD13Q1 (V006) EVI data from the official website of the United States Geological Survey (USGS) (https://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.006/). The data is synthesized from 16 days of data. There are 23 episodes per year, and the spatial resolution is 250m×250m.This paper takes Hunan Province as the experimental area. Hunan Province covers the image data of H27V05, H27V06, and H28V06. During the 21 years, there were a total of 483 episodes with 1449 images. The MOD13Q1 data has a total of 12 bands. First, MRT (Modis Reprojection Tool) software is used to extract EVI data and quality file, which is applied to image mosaic, format conversion, projection transformation, resampling and other preprocessing. Then, the preprocessed EVI file and quality file in the administrative division vector data of Hunan Province are extracted by mask to obtain the EVI spatio-temporal data and quality file spatio-temporal data of Hunan Province from 2001 to 2021.

2.2 Vegetation remote sensing time series data reconstruction

Clouds, fog, ice and snow, and other factors can bring noise to the MODIS 13Q1 data, resulting in unreliable or null values in the EVI time series data. The EVI pixels with unreliable quality are replaced by the pixel values reconstructed from EVI time series data Savitzky-Golay filter sequence at the pixel position. Invalid fields are temporarily set to null. As for the null value, the method of replacing the time-series data with the average value of the image in this period for many years or the average value of the 3*3 spatial neighborhood in the same period is used for filling. Finally, a total of 483 EVI image time series data were reconstructed and outliers removed. Figure 1 shows the processed EVI images of 2001 and 2021 on the 209th day (July 28).

Figure 1. EVI images of Hunan Province on July 28, 2001 and 2021 after processing

2.3 EVI time series data acquisition of remote sensing images

According to the basic principle of spatio-temporal series prediction, the EVI image data of the study area are transformed in spatial dimension and structure. The image time series data representing the two-dimensional continuous space vegetation change is converted into a one-dimensional area center point vegetation time series data array representing the subspace vegetation change. Among them, the formed one-dimensional point time series data meets the input requirements of the LSTM prediction model.

Figure 2. Distribution of sample points and test points of EVI time series data of vegetation change prediction in Hunan Province

There are 3.9363 million pixels point time series within the effective range of the boundary of the study area. Due to equipment conditions limitations, in this paper, the systematic regular sampling method is adopted for thinning sampling. A total of 788 sampling points regularly distributed in the study area are obtained. The EVI time series data of the corresponding positions in the 483 EVI images are assigned to the sample points, with each sample time series data length of 483. A total of 788 sample point vegetation time series data sets is used as the data source of vegetation change prediction. On the basis of 788 sample points, 8 points are randomly selected to test the performance of different time scale sequences in prediction after discrete wavelet decomposition, and to further evaluate the performance of the method in the overall sample. The distribution of EVI time series data sample points and test points are shown in Figure 2. The label of the test point is the number of the point in the EVI time series sample point set.

3. Research Methods

3.1 Principle of discrete wavelet decomposition

The wavelet transform (WT) is a method of multi-resolution analysis in time and frequency domains. It has been widely used in the fields of meteorology, hydrology and geophysics [30-33]. The basic idea of wavelet transform is the decomposition of a signal into different time scales with a set of basis functions. The set of basis functions {ψa,b(t)} can be generated by translating and scaling the so-called mother wavelet function ψ(t), according to the Equation 1:

$\begin{gathered}\psi_{a, b}(t)=\frac{1}{\sqrt{a}} \psi\left(\frac{t-b}{a}\right), \\ a>0, \quad-\infty<b<\infty\end{gathered}$          (1)

where, a is the scale parameter which adjusts the dilation of the wavelet and b determines the location of the wavelet. For a time-series f(t), the Continuous Wavelet Transform (CWT) of f(t) can be defined as Equation 2:

$W_\psi(a, b)=\frac{1}{\sqrt{a}} \int_{-\infty}^{+\infty} \bar{\psi}\left(\frac{t-b}{a}\right) f(t) d t$          (2)

where, $\bar{\psi}$ is the mother wavelet complex conjugate. For each scale a, the result of the WT is a set of coefficients Wψ(a,b), associated with different locations b.

In real applications, the continuous wavelet is usually discretized into discrete signal for processing. In Formula (1), we take a=2j, b=k2j, where k is the location index and j is referred to as the decomposition level. Thus, a discretely scaled and translated wavelet basis can be expressed as Equation 3:

$\psi_{j, k}(t)=2^{-j / 2} \psi\left(2^{-j} t-k\right)$           (3)

and the discrete wavelet transform (DWT) of f(t) can be written as Equation 4:

$\begin{aligned} & W_\psi(j, k)=2^{-j / 2} \int_{-\infty}^{+\infty} \bar{\psi}\left(2^{-j} t-k\right) f(t) d t, \\ & j=0,1, \ldots, k \in Z\end{aligned}$           (4)

The characteristics of the origin time series f(t) at the decomposition j and time location k can be represented by Wψ(j,k). At a lower level j, the smaller and finer component of the signal is accessed.

Multi-resolution analysis (MRA) based on DWT decomposes a signal with different frequencies into a certain number of sub-signals at different time scales by successively translating and convolving the elements of a sequence high-pass filter and low-pass filter associated with the mother wavelet. These filters preserve the details of the signal (D) and the approximate components (A) of the sub sequence. In the first level of the decomposition, f(t)=A1+D1, the signal has a low-pass filtered component A1, and a high-pass filtered component D1. In order to obtain a decomposition on a coarser scale, A1=A2+D2, the same procedure is performed on A1. Thus, we have Aj-1=Aj+Dj for j=2,...,J. Then an original time series, which is decomposed into J levels by DWT, can be expressed as follows:

$f(t)=\mathrm{D}_1(\mathrm{t})+\mathrm{D}_2(\mathrm{t})+\ldots+\mathrm{D}_J(\mathrm{t})+\mathrm{A}_J(\mathrm{t})$           (5)

where, J is the highest decomposition level considered. The detail component Dj at a particular decomposition level j is given by Equation 6.

$D_j(t)=\sum_{k=-\infty}^{\infty} W_\psi(j, k) \psi_{j, k}(t)$          (6)

The approximation component Aj of the signal on the scale 2j is written as Equation 7:

$A_j(t)=\sum_{k=-\infty}^{\infty} V_{\varphi}(j, k) \varphi_{j, k}(t)$           (7)

where, ϕj,k(t) is the scaling function and Vϕ(j,k) are the scaling coefficients computed by ϕj,k(t). Suppose the sum of detail changes after J level discrete wavelet decomposition is VJ(t), then VJ(t) can be calculated by the Equation 8:

${{V}_{J}}(t)=\sum\limits_{j=1}^{J}{{{D}_{j}}}$           (8)

According to Equation 5 and Equation 8, the original time series can be decomposed into approximate components and the sum of detailed components of different time scales, so f(t) can be expressed in the following form in Equation 9:

$f(t)={{A}_{1}}(t)+{{V}_{1}}(t)={{A}_{2}}(t)+{{V}_{2}}(t)=...={{A}_{J}}(t)+{{V}_{J}}(t)$            (9)

Therefore, the original time series can be regarded as the sequence combination of trend and detailed change components on different time scales.

3.2 Long short term memory networks prediction model

Long short term memory networks (LSTM) is a recurrent neural network (RNN) improved by gating proposed by Hochreiter and Schmidhuber [34] in 1997, which is used to solve the problem of gradient disappearance or gradient explosion. The model introduces gated recurrent unit (GRU) and long short term memory unit on the basis of RNN. This model captures deep connections better and improves the gradient disappearance problem by changing the hidden layer of RNN [28].

Similar to RNN, LSTM recurrent network structure is composed of multiple layers of cells, including input layer, output layer and hidden layer. The difference between LSTM and RNN lies in the design of the hidden layer. After the internal processing of the hidden layer cell, the hidden layer cell of LSTM has two transmission states, which are the cell state representing the current neural cell to store long-term information and the hidden state representing the current neural cell to store short-term information. The node cell structure of LSTM is shown in Figure 3. The node cell of the LSTM loop network is generally composed of input gate, forget gate, output gate and cell state.

Figure 3. LSTM neural network cell structure diagram

The control functions of each gate and cell state in Figure 3 are as follows:

Forget gate ${{f}_{t}}=\sigma \left( {{W}_{f}}{{x}_{t}}+{{U}_{f}}{{h}_{t-1}}+{{b}_{f}} \right)$           (10)

Input gate ${{i}_{t}}=\sigma \left( {{W}_{i}}{{x}_{t}}+{{U}_{i}}{{h}_{t-1}}+{{b}_{i}} \right)$          (11)

Newly added state data ${{\tilde{c}}_{t}}=\tanh \left( {{W}_{c}}{{x}_{t}}+{{U}_{c}}{{h}_{t-1}}+{{b}_{c}} \right)$           (12)

Cell state ${{c}_{t}}={{f}_{t}}{{c}_{t-1}}+{{i}_{t}}{{\tilde{c}}_{t\text{ }}}$           (13)

Output gate ${{o}_{t}}=\sigma \left( {{W}_{o}}{{x}_{t}}+{{U}_{o}}{{h}_{t-1}}+{{b}_{o}} \right)$          (14)

Hidden state ${{h}_{t}}={{o}_{t}}\tanh ({{c}_{t}})$          (15)

In Equations (10)~(14), xt is the input data at the current moment and ht-1 is the hidden state output value of the cell state at the previous moment. These two data are simultaneously input to the forget gate, input gate, output gate and cell state functions in the LSTM node cell to jointly determine the output values of these functions. In Equations (10)~(12) and (14), W, U and b represent the input weight, loop weight, and bias of each gate and cell state, which are different in each layer. They are combined with the input variables xt and ht-1 linear transformation and then mapped by the activation function sigmoid (σ) or tanh cell to obtain the corresponding weight value. As for the activation function, the σ function and the tanh function are as follows:

$sigmoid(x)=\frac{1}{1+{{e}^{-x}}}$           (16)

$\tanh (x)=\frac{{{e}^{x}}-{{e}^{-x}}}{{{e}^{x}}+{{e}^{-x}}}$            (17)

In order to control the degree of forgetting, memory, output, hidden state and newly added state data, the Sigmoid function ensures that the weight value of the output of the gating mechanism is between [0,1]. The tanh function is a hyperbolic tangent function, which ensures that the element values of the state are between [-1,1]. When the output gate is approximately 1, the information contained in the cell state will be output. When the output gate is approximately 0, the information contained in the cell state will be preserved.

4. Experiment

4.1 Discrete wavelet multi-scale decomposition of remote sensing image EVI time series data

According to Martinez and Gilabert [3] and Long et al. [5], this paper adopts discrete Meyer wavelet (dmey) for decomposition, which has biorthogonal and compact support, and is usually used for fast discrete wavelet transformation. In order to select the most suitable scale for the study of the interannual component and intra-annual components of vegetation change, the period or half-period is calculated at different scales of the discrete Meyer wavelet MRA decomposition. The calculation method is shown in Eq. (18):

$p=\frac{a\Delta t}{{{V}_{c}}}$          (18)

In Eq. (18), a is the time scale, Δt is the sampling period, Vc is the center frequency of the wavelet base, and p is the period under the corresponding scale of MRA.

According to the discrete Meyer wavelet base center frequency Vc=0.6634 given by MATLAB R2018b and the vegetation index sampling period Δt=16 days, the half-period and period of each scale are calculated. As shown in Table 1, when the time scale a=25, the corresponding half-cycle p/2=386 days. Since the optimal time scale of wavelet transform is consistent with the half-cycle, the approximation component for level 5th, A5 provides information about the inter-annual variability over the study period. This approximation component signals how much change has occurred over time scales of 386 days or more, enabling identification of the long-term variation in the study period. Similarly, A3 and A1 represent approximate components of seasonal scale and monthly scale respectively, and the variation trend of EVI at these two scales can be seen. Conversely, the detail components D1~D5, with semi-periods ranging from 24 to 386 days, are used to analyze the intra-annual vegetation dynamics, including monthly, quarterly, half-year scale changes. D1, D2, D3, D4, D5 represent the detailed changes of vegetation dynamics on the scale of 1-24 days, 24-48 days, 48-96.5 days, 96.5-193 days, and 193-386 days, respectively.

The test point 85 (latitude of 112.3809° and longitude of 29.3044°) located in the savanna vegetation type in the Dongting Lake area is selected for discrete wavelet decomposition, and the decomposition results are shown in Figure 4. According to the A2 composition, it can be seen that there are multiple peaks in this pixel year. Among the A3 trend components representing the seasonal scale, there is an obvious EVI peak every year. According to D3, a more detailed change in the “M-shaped” seasonal scale can be seen. As the time scale grows, the finer components of the signal representing intra-annual variation are preserved in the detail component sequence (right column of Figure 4), and the smoother part is given by an approximate component sequence associated with the low-frequency components (left column of Figure 4). It can be seen that EVI data is subjected to wavelet multi-scale decomposition and reconstruction to obtain approximate component information and detailed component information at different time scales. MRA is used for temporal multi-scale decomposition of EVI image time series data. MATLAB is used to read 483 EVI remote sensing data from 2001 to 2021 into a three-dimensional array. Through the discrete Meyer wavelet transform, the data in the study area were subjected to 5-level discrete wavelet decomposition and reconstruction pixel by pixel, so as to obtain the A5, D5, D4, D3, D2, D1 subsequence and complete the decomposition.

Table 1. Period (p) and half-period (p/2) of discrete Meyer DWT under different decomposition levels of MRA with a sampling period of 16 days

level (j)

scale (a)

p(period days)

p/2(half-period days)

1

2

48

24

2

4

96.5

48

3

8

193

96.5

4

16

386

193

5

32

772

386

Figure 4. Multi-resolution analysis of the EVI time series at test point 85 (The original time series and 5 approximate components are shown on the left, and the 5 detail components are shown on the right)

4.2 Generation of multiscale vegetation EVI decomposition sequence group

On the basis of the subsequence obtained by discrete wavelet 5-level decomposition in the above section, the wavelet decomposition principle is used to calculate the approximate component A4, A3, A2, A1 of the other four levels of decomposition and the sum of the detailed components of each level of decomposition, namely V5, V4, V3, V2, V1. In this section, according to the law of vegetation change, the monthly scale (A1 and V1), seasonal scale (A3 and V3), and annual scale (A5 and V5) sequence group images are obtained, and the sampling points assigned to the system by A1 and V1, A3 and V3, and A5 and V5 time series data are extracted. Among them, the sequence group data of 8 test points are extracted. Figure 5 shows the original EVI time series data at test point 85 and its decomposed A1 and V1, A3 and V3, A5 and V5 component sequence data. After decomposing, the sum of the approximate component A and the detail component V is equal to the value of the EVI original time series data. In other words, the time series conforms to the relation EVI=A1+V1=A3+V3=A5+V5, and the wavelet multi-scale decomposition is correct.

The data statistics data of EVI, A1, A3, A5, V1, V3, and V5 sequence data at this point are shown in Table 2. It can be seen from Table 2 that the maximum value of each scale trend component becomes smaller, the minimum value becomes larger, and the value range and standard deviation become smaller with the increase of monthly scale, seasonal scale and annual scale. As for the detail components of the corresponding monthly scale, seasonal scale, annual scale, the maximum value increases and the minimum value decreases. Furthermore, more detail components are separated out. As the decomposition scale increases, the range and standard deviation of the sum of detail components increase. The decomposed approximate components and detailed component sequences behave exactly the opposite way, but the mean value of the EVI raw series and its trend component are the same, and the mean value of the detail component tends to zero, and their performance in EVI prediction remains to be explored. Next, 8 test points are used for multi-time scale vegetation time series data prediction experiments.

Figure 5. Schematic diagram of the original EVI data and its wavelet decomposition sequences (test point 85)

Table 2. Statistical table of EVI time series and its decomposition sequences data (test point 85)

 

EVI

A1

A3

A5

V1

V3

V5

max

0.5563

0.5632

0.4092

0.2784

0.0653

0.1974

0.3068

min

0.0155

0.0327

0.0657

0.2016

-0.0721

-0.1867

-0.2349

range

0.5408

0.5305

0.3435

0.0769

0.1374

0.3841

0.5417

mean

0.2486

0.2486

0.2487

0.2486

0.0000

-0.0001

0.0000

Stdev

0.1203

0.1191

0.0945

0.0163

0.0169

0.0744

0.1191

4.3 Multi-scale EVI time series group LSTM modeling

The LSTM prediction comparison of monthly scale, seasonal scale, annual scale sequence group data and original time series data is carried out to verify whether the decomposition sequence prediction is better than the original time series prediction, and to find the optimal decomposition scale. First, the time series groups A1 and V1,A3 and V3, A5 and V5, and the original EVI time series are tested and evaluated by LSTM modeling at the 8 test points. Then, the optimal decomposition scale sequence is used for 788 sample point predictions. The overall prediction accuracy is evaluated by comparing it with the original time series prediction accuracy.

Figure 6. LSTM modeling prediction results of monthly scale, seasonal scale, annual scale sequence group and original EVI time series in test point 85

In order to compare the prediction effect of different scale sequence group data and the original EVI time series data change, LSTM, which is excellent in spatiotemporal series prediction, is used for modeling prediction. The construction of the LSTM model adopts the Pytorch deep learning framework. 8 test points were tested multiple times for LSTM model building. After comparative analysis, the main parameters are set as follows. The hidden layer dimension (hidden_size) is set to 128, the number of trainings (num_epochs) is set to 10000, and the batch size (batch_size) is set to 64. As for the learning rate (learning_rate), it is set to 0.0004 for A1, A3, A5 approximation component prediction and 0.0002 for V1, V3, V5 detail component and EVI original time series prediction. The sliding window length is set to 69 for approximate component prediction and 23 for detail component and EVI original time series prediction. The model optimizer chooses Adam, and the dropout is set to 0.7. In the time series data, a total of 368 time point data from 2001 to 2016 are set as the training set Xtrain, and a total of 115 time point data from 2017 to 2021 are set as the test set Xtest. The LSTM with the parameter settings mentioned above is used to model and predict the time series of A1, A3, A5, V1, V3, V5, and EVI of the 8 test points. Figure 6 shows the decomposition sequence of each time scale in the test point 85 and the prediction results of the original time series.

4.4 Analysis of prediction results of multi-scale EVI sequence group

Table 3 counts the loss rate predicted by LSTM modeling of A1, A3, A5, V1, V3, V5, and the original EVI time series at 8 test points. Among them, Loss_EVI, Loss_A1, Loss_A3, Loss_A5, Loss_V1, Loss_V3, and Loss_V5 are the predicted loss rates of the original EVI time series and sequences A1, A3, A5, V1, V3, V5, respectively.

Table 3. Test point decomposition sequence and original time series LSTM prediction loss rate statistics (retain four decimal places)

Point

ID

Loss_

EVI

Loss_

A1

Loss_

A3

Loss_

A5

Loss_

V1

Loss_

V3

Loss_

V5

85

0.0019

0.0008

0.0001

0.0000

0.0000

0.0012

0.0026

95

0.0113

0.0041

0.0002

0.0000

0.0004

0.0052

0.0110

310

0.0055

0.0010

0.0001

0.0000

0.0001

0.0027

0.0054

372

0.0089

0.0027

0.0002

0.0000

0.0002

0.0042

0.0089

407

0.0058

0.0024

0.0002

0.0000

0.0001

0.0031

0.0057

428

0.0058

0.0027

0.0003

0.0001

0.0001

0.0025

0.0053

605

0.0027

0.0015

0.0001

0.0000

0.0001

0.0012

0.0025

707

0.0048

0.0014

0.0001

0.0000

0.0002

0.0028

0.0041

Mean

0.0059

0.0021

0.0002

0.0000

0.0002

0.0029

0.0057

According to the loss rate statistics (Table 3), the predicted loss rate is mostly within 10‰, and the maximum value is controlled at about 11.3‰. Thus, the accuracy meets the requirements. The larger the time scale, the more approximate components the detail V contains, the greater the prediction loss rate of the detail component, and the greater the error of the prediction result. At the same time, the less detailed component the approximate component A contains, the smaller the error of the predicted result and the more accurate the prediction. The loss rate mean is found. The loss rate of V5 is close to EVI, and the loss rates of other sequences are much smaller than the EVI loss rate. By summing the loss rate of A1 and V1, A3 and V3, A5 and V5 and comparing it with the loss rate of the original time series of EVI, it can be known that Loss_A1+Loss_V1<Loss_A3+Loss_V3<Loss_A5+Loss_V5<Loss_EVI. The sum of loss rates for all scaled sequence predictions is smaller than the original time series prediction.

Figure 7. Monthly scale, seasonal scale, annual scale prediction error and original EVI time series error of test point 85

Table 4. Comparison of prediction accuracy between the monthly scale, seasonal scale, and annual scale sequence groups at the test points and the EVI original time series

pointID

Series

R2

MAE

RMSE

TIC

85

A1+V1

0.9430*

0.0222*

0.0286*

0.0506*

A3+V3

0.9108

0.0284

0.0358

0.0635

A5+V5

0.8178

0.0387

0.0512

0.0895

EVI

0.8678

0.0328

0.0436

0.0781

 

 

 

 

 

 

95

A1+V1

0.8392*

0.0520*

0.0664*

0.0749*

A3+V3

0.8030

0.0571

0.0735

0.0836

A5+V5

0.5984

0.0847

0.1049

0.1190

EVI

0.5881

0.0845

0.1063

0.1214

 

 

 

 

 

 

310

A1+V1

0.9288*

0.0258*

0.0334*

0.0457* 

A3+V3

0.8192

0.0415

0.0532

0.0737

A5+V5

0.6548

0.0557

0.0735

0.1014

EVI

0.6464

0.0560

0.0744

0.1032

 

 

 

 

 

 

372

A1+V1

0.8689*

0.0420*

0.0544*

0.0611*

A3+V3

0.8032

0.0524

0.0666

0.0758

A5+V5

0.6065

0.0740

0.0942

0.1074

EVI

0.6053

0.0740

0.0943

0.1085

 

 

 

 

 

 

407

A1+V1

0.8723*

0.0389*

0.0497*

0.0657*

A3+V3

0.8315

0.0455

0.0571

0.0756

A5+V5

0.7092

0.0591

0.0750

0.0978

EVI

0.6979

0.0601

0.0765

0.1013

 

 

 

 

 

 

428

A1+V1

0.8506

0.0430

0.0539

0.0722

A3+V3

0.8540*

0.0428*

0.0533*

0.0720*

A5+V5

0.7314

0.0579

0.0723

0.0970

EVI

0.7021

0.0609

0.0762

0.1054

 

 

 

 

 

 

605

A1+V1

0.9137

0.0313

0.0388

0.0617

A3+V3

0.9211*

0.0294*

0.0371*

0.0586*

A5+V5

0.8550

0.0386

0.0503

0.0805

EVI

0.8436

0.0416

0.0522

0.0840

 

 

 

 

 

 

707

A1+V1

0.8295*

0.0295*

0.0373*

0.0494*

A3+V3

0.6480

0.0408

0.0536

0.0710

A5+V5

0.4893

0.0486

0.0646

0.0859

EVI

0.4100

0.0528

0.0694

0.0926

Note: *Bold marks are the best indicator results for each point. The oblique bold mark means that the prediction index result is less accurate than the original time series prediction result.

The LSTM prediction results of the 2017-2021 test set data of A1, A3, A5, V1, V3, V5, and EVI time series of the 8 test points were counted separately. The monthly scale prediction is obtained by adding the prediction of A1 and V1. The seasonal scale prediction is obtained by adding the A3 and V3 prediction. The annual scale prediction is obtained by adding the A5 and V5 predictions. The LSTM prediction results of the monthly scale, seasonal scale, and annual scale sequences and original EVI time series of the 8 test points, and the errors between these time series prediction results and the original EVI data are obtained. There are 4 error values in total, which are error 1, error 3, error 5, and error raw. The four errors of test point 85 are shown in Figure 7. Among them, the monthly scale decomposition sequence prediction result of point 85 is the best. The error order is monthly scale<seasonal scale<annual scale<original time series.

The prediction accuracy evaluation indicators R2, MAE, RMSE, and TIC of the monthly scale, seasonal scale, annual scale, and EVI original time series of the eight test points are calculated and the results are summarized in Table 4. R2, MAE, RMSE, and TIC are coefficient of determination, mean absolute error, root mean square error, and Hill inequality coefficient, respectively.

In Table 4, the performance of the four accuracy indicators on the 8 test points are respectively counted. R2 represents the correlation with the original time series. The closer R2 is to 1, the higher the accuracy. As for the three indicators of MAE, RMSE, and TIC that represent the error, the closer the value is to 0, the higher the accuracy. Among the 8 test points, there are 6 test points (test point 85, test point 95, test point 310, test point 372, test point 407 and test point 707) with the largest R2 value and the smallest MAE, RMSE, and TIC values in the monthly scale sequence group prediction. The prediction accuracy of monthly scale is higher than that of seasonal scale, annual scale, and undecomposed original time series. In summary, the order of prediction accuracy is monthly scale > seasonal scale > annual scale > original time series. The monthly scale has the best performance. Another two test points (test point 428 and test point 605) have the largest R2 value and the smallest MAE, RMSE, and TIC values in the prediction of the seasonal scale sequence group. The prediction accuracy of seasonal scale is better than that of monthly scale, annual scale and the original time series. In summary, the order of prediction accuracy is seasonal scale > monthly scale > annual scale > original time series.

With the increase of monthly scale, seasonal scale, and annual scale, the R2 values of most test points are getting lower and lower, and the MAE, RMSE, and TIC values are getting higher and higher. In general, the prediction accuracy of the decomposed sequence group is higher than that of the original time series. For individual test points, such as test point 85, the prediction results on the annual scale decomposition sequence are worse than the direct prediction accuracy of the original time series, which may be related to the fixed LSTM parameters used in the experiment. This parameter may not be well suited for predicting at the annual scale of the time series data at that point. At the same time, the annual scale decomposition sequence prediction results of test points 95 and 372 are greater than the original time series prediction results in terms of MAE indicators. However, the two sequence forecast accuracy metrics are close in numerical value. Therefore, it can be concluded that the prediction results of the annual scale sequence group are unstable. In contrast, both seasonal scale and monthly scale decomposition sequence prediction results are better than the original time series, with relatively stable performance.

Table 5 further performs mean statistics on the accuracy evaluation indicators of these 8 test points. The monthly scale prediction performance index represented by the A1+V1 sequence is average and the overall performance is the best. The mean values of the precision evaluation indexes of A3+V3 and A5+V5 are better than the original EVI time series. The order of accuracy is monthly scale>seasonal scale>annual scale>original EVI time series. This also verifies that the sequence after discrete wavelet decomposition can facilitate the learning of approximate and detailed component features and improve prediction performance.

Table 5. Comparison of the mean value of the prediction accuracy index between the test point multi-scale sequence group and the EVI original time series

Index Item

A1+V1

A3+V3

A5+V5

EVI

R2

0.8808

0.8238

0.6828

0.6701

MAE

0.0356

0.0422

0.0572

0.0578

RMSE

0.0453

0.0538

0.0732

0.0741

TIC

0.0601

0.0717

0.0973

0.0993

Note: The numbers marked in black and bold are the best index results for each accuracy index item.

According to the results of the test point prediction experiment, the A1 and V1 decomposition sequence groups with the best prediction performance and the A3 and V3 decomposition sequence groups with the better prediction performance are selected to perform LSTM modeling and prediction on 788 sample points, and then compared with the prediction accuracy of the original EVI time series. The mean value of the accuracy evaluation indicators of all points is calculated to Table 6. It can be seen from Table 6 that the accuracy indicators of the monthly scale decomposition sequences A1 and V1 of 788 sample points are better than the seasonal scale decomposition sequences and original time series modeling predictions.

Table 6. Comparison of the average prediction accuracy between the monthly and seasonal scale sequence group of all sample points and the original time series of EVI

Index Item

A1+V1

A3+V3

EVI

R2

0.8852

0.7761

0.6741

MAE

0.0346

0.0477

0.0568

RMSE

0.0438

0.0607

0.0733

TIC

0.0616

0.0853

0.1046

5. Conclusions

Through discrete wavelet transformation, the MODIS EVI data series with a sampling interval of 16 days in Hunan Province from 2001 to 2021 are decomposed into 5 levels pixel by pixel. The monthly scale, seasonal scale, and annual scale vegetation approximate components and detailed component sequences of each pixel are obtained. Based on the decomposed multi-scale EVI time series data, the LSTM modeling prediction and accuracy evaluation comparative analysis of EVI time series data monthly scale, seasonal scale, and annual scale discrete wavelet decomposition sequences are completed.

For the collected 16d time resolution EVI data, the discrete wavelet decomposition method can effectively decompose the vegetation change sequence of approximate monthly scale, seasonal scale, and annual scale. Each scale sequence group of wavelet decomposition strips off trend features and detail features. The decomposed sequence group can still reconstruct the original time series and retain all the information of the original time series. Therefore, the prediction of the sequence group can reflect the prediction of the original time series.

In this paper, the monthly, seasonal and annual scale sequences depict the characteristics of vegetation change at different scales respectively. In the discrete wavelet decomposition, each high-level detailed components contains high-frequency information of each level and part of low frequency information of the lower level, and each high-level approximate components removes more high frequency information of the original data. Therefore, the higher the number of decomposition level, the greater the prediction loss of high-frequency detail information, and the lower the loss of low-frequency trend information. The different sequence groups have different predictive ability and efficiency. In general, the vegetation change prediction effect of the sequence group data decomposed by EVI monthly scale, seasonal scale, and annual scale is better than that of the original EVI time series. Discrete wavelet decomposition helps to improve the accuracy of vegetation change prediction. According to the experimental results, the ranking of the vegetation change prediction performance of the sequence group decomposed by monthly scale, seasonal scale, and annual scale is: monthly scale > seasonal scale > annual scale > original EVI time series.

In the prediction experiment, the original EVI time series, monthly scale decomposition sequences (A1 and V1) and seasonal scale decomposition sequences (A3 and V3) of 788 sample points are respectively subjected to LSTM modeling and vegetation time series data prediction point by point. Although this method takes into account the spatial heterogeneity of vegetation change, it takes a long time. The time predicted by the monthly scale series was almost twice as long as that predicted by the original EVI series. Therefore, the next step is to consider building a model that takes into account both the temporal and spatial heterogeneity of sample time series data and the time complexity of prediction modeling for large-scale vegetation change prediction. For example, zoning or classification modeling is used to reduce the number of models so as to reduce the total time complexity.

Acknowledgments

This work is supported by the General Project of Education Department of Hunan Province (Grants No.: 20C0373, 22C0500), the Office of leading group of Education Science Planning of Guangdong Province (Grant No.: 2019GXJK105), the Education Department of Guangdong Province (Grant No.: 2020WTSCX126), and the Key Laboratory of Remote Sensing Monitoring of Dongting Lake (Grants No.: DTH Key Lab.2022-09, DTH Key Lab.2022-06).

  References

[1] Chi, J., Lee, H., Hong, S.G., Kim, H.C. (2021). Spectral characteristics of the antarctic vegetation: A case study of barton peninsula. Remote Sensing, 13(13): 2470. https://doi.org/10.3390/rs13132470

[2] Robert, N.C. (1967). Remote sensing as a means of determining ecological conditions. BioScience, 17: 444-449. https://doi.org/10.2307/1293814

[3] Martinez, B., Gilabert, M.A. (2009). Vegetation dynamics from NDVI time series analysis using the wavelet transform. Remote Sensing of Environment, 113: 1823-1842. https://doi.org/10.1016/j.rse.2009.04.016

[4] Jin, Y.S., Jin, K., Wang, F., Liu, C.X., Qin, P., Zong, Q.L., Liu, P.R., Chen, M.L. (2022). Imapacts of climate and human activities on NDVI change in eastern coastal areas of China. Environmental Science, 1-18. https://doi.org/10.13227/j.hjkx.202207039

[5] Long, Y.H; Qin, J.X., He, X.G., Yang, Z. (2015). Wavelet multi-resolution analysis of vegetation dynamic change in Dongting Lake Basin. Acta Geographica Sinica, 70: 1491-1502. https://doi.org/10.11821/dlxb201509011

[6] Hao, H.G., Li, Y.Y., Zhang, H.Y., Zhai, R.X., Liu, H.Y. (2019). Spatio-temporal variations of vegetation and its determinants in the National Key Ecological Function Area on Loess Plateau between 2000 and 2015. Ecology and Evolution, 9: 5810-5820. https://doi.org/10.1002/ece3.5165

[7] Ding, Q., Liang, Y.J., Xu, N.H., Nai, J.M., Hu, X.M. (2022). Inversion of soil moisture by BDS dual-frequency multi-star combination under low vegetation cover environment. Progress in Geophysics (in Chinese), 37: 2242-2250. https://doi.org/10.6038/pg2022FF0525

[8] Guijarro, M., Riomoros, I., Pajares, G., Zitinski, P. (2015). Discrete wavelets transform for improving greenness image segmentation in agricultural images. Computers and Electronics in Agriculture, 118: 396-407. https://doi.org/10.1016/j.compag.2015.09.011

[9] Pang, J.F., Liu, B., Zhang, B., Zhang, P.P., Wang, B., Zeng, F.J. (2019). Study on the combined model of forecasting the days of sand-dust weather based on wavelet decomposition--taking the time series of dust weather in the trasitional zone of qira desert oasis during 2008-2016 as an example. Metor Mon, 45: 651-658. https://doi.org/10.7519/j.issn.1000-0526.2019.05.006

[10] Yan, J.H. (2017). Stock index portfolio prediction based on discrete wavelet decomposition and support vector machine. Statistics & Decision, 471: 168-171. https://doi.org/10.13546/j.cnki.tjyjc.2017.03.041

[11] Jin, B.B., Hu, Y. (2022). Spatial-temporal wavelet analysis video prediction based on differential attention mechanism. Journal of Computer-Aided Design & Computer Graphics, 34: 180-188. https://doi.org/10.3724/SP.J.1089.2022.18895.

[12] Paul, R.K., Garai, S. (2022). Wavelets based artificial neural network technique for forecasting agricultural prices. Journal of the India Society for Probability and Statistics, 23: 47-61. https://doi.org/10.1007/s41096-022-00128-3

[13] Rhif, M., Ben Abbes, A., Martinez, B., Farah, I.R. (2021). An improved trend vegetation analysis for non-stationary NDVI time series based on wavelet transform. Environmental Science and Pollution Research, 28: 46603-46613. https://doi.org/10.1007/s11356-020-10867-0

[14] Rhif, M., Ben Abbes, A., Martinez, B., Farah, I.R., Gilabert, M.A. (2022). Optimal selection of wavelet transform parameters for spatio-temporal analysis based on non-stationary NDVI MODIS time series in Mediterranean region. ISPRS Journal of Photogrammetry and Remote Sensing, 193: 216-233. https://doi.org/10.1016/j.isprsjprs.2022.09.007

[15] Shi, S., Li, W., Ding, Y.S., Lin, X.P., Zhai, Y.C. (2023). Spatiotemporal evolution and influencing factors of vegetation in Northeast China. China Environmental Science, 43: 276-289. https://doi.org/10.19674/j.cnki.issn1000-6923.2023.0008

[16] Guo, A.P., Jiang, A.J., Lin, J., Li, X.X. (2020). Data mining algorithms for bridge health monitoring: Kohonen clustering and LSTM prediction approaches. The Journal of Supercomputing, 76: 932-947. https://doi.org/10.1007/s11227-019-03045-8

[17] Liu, J.W., Wang, Y.F., Luo, X.L. (2021). Research and Development on Deep Memory Network. Chinese Journal of Computers, 44: 1549-1589. https://doi.org/10.11897/SP.J.1016.2021.01549

[18] Yang, Q., Qin, L., Gao, P., Zhang, R.B. (2021). Prediction of annual precipitation in the Northern Slope Economic Belt of Tianshan Mountains based on a EEMD-LSTM model. Arid Zone Research, 38: 1235-1243. https://doi.org/10.13866/j.azr.2021. 05. 05

[19] Bhandari, H.N., Rimal, B., Pokhrel, N.R., Rimal, R., Dahal, K.R., Khatri Rajendra, K.C. (2022). Predicting stock market index using LSTM. Machine Learning with Applications, 9: 100320. https://doi.org/10.1016/J.MLWA.2022.100320

[20] Lin, H.L., Sun, Q.B., Chen, S.Q. (2020). Reducing Exchange Rate Risks in International Trade: A Hybrid Forecasting Approach of CEEMDAN and Multilayer LSTM. Sustainability, 12: 2451. https://doi.org/10.3390/su12062451

[21] Zhao, Z., Chen, W.H., Wu, X.M., Chen, P.C.Y., Liu, J.M. (2017). LSTM network: A deep learning approach for short-term traffic forecast. IET Intelligent Transport Systems, 11: 68-75. https://doi.org/10.1049/iet-its.2016.0208

[22] Wang, J.Q., Du, Y., Wang, J. (2020). LSTM based long-term energy consumption prediction with periodicity. Energy, 197: 117197. https://doi.org/10.1016/j.energy.2020.117197

[23] Murugesan, R., Mishra, E., Krishnan, A.H. (2022). Forecasting agricultural commodities prices using deep learning-based models: basic LSTM, bi-LSTM, stacked LSTM, CNN LSTM, and convolutional LSTM. International Journal of Sustainable Agricultural Management and Informatics, 8: 242-277. https://doi.org/10.1504/IJSAMI.2022.10048228

[24] Shen, H.J., Luo, Y., Zhao, Z.C., Wang, H.J. (2020). Prediction of summer precipitation in China based on LSTM network. Climate Change Research, 16: 263-275. https://doi.org/10.12006/j.issn.1673-1719.2019.067

[25] Harleen, K., Afshar, A.M., Saleha, M., Bhavya, A., Ritu, C., Muhammad, A.R., Ozgur, K. (2021). Predicting water availability in water bodies under the influence of precipitation and water management actions using VAR/VECM/LSTM. Climate, 9: 144. https://doi.org/10.3390/cli9090144

[26] Bao, R.Y., He, Z.H., Zhang, Z.Y. (2022). Application of lightning spatio-temporal localization method based on deep LSTM and interpolation. Measurement, 189: 110549. https://doi.org/10.1016/j.measurement.2021.110549

[27] Drewil, G.I., Al-Bahadili, R.J. (2022). Air pollution prediction using LSTM deep learning and metaheuristics algorithms. Measurement: Sensors, 24: 100546. https://doi.org/10.1016/J.MEASEN.2022.100546

[28] Vasilakos, C., Tsekouras, G.E., Kavroudakis, D. (2022). LSTM-based prediction of mediterranean vegetation dynamics using NDVI time-series data. Land, 11(6): 923. https://doi.org/10.3390/land11060923

[29] Reddy, D.S., Prasad, P.R.C. (2018). Prediction of vegetation dynamics using NDVI time series data and LSTM. Modeling Earth Systems and Environment, 4: 409-419. https://doi.org/10.1007/s40808-018-0431-3

[30] Sang, Y.F., Wang, Z.G., Liu, C.M. (2013). Applications of wavelet analysis to hydrology: Status and prospects. Progress in Geography, 32(9): 1413-1422. https://doi.org/10.11820/dlkxjz.2013.09.011

[31] Wang, Z.T., Tao, Y.W., Zhang, A.L., Wang, Y.K., Wang, D. (2022). Analysis of time series characteristics of typical hydrometeorological elements in Diaojiang River basin mining area in Guangxi. Journal of China Hydrology, 42: 83-87. https://doi.org/10.19797/j.cnki.1000-0852.20210306

[32] Tian, J.H., Li, Z.Q., Ke, H.C., Liang, L., Tian, M.H., Gao, J.F. (2023). Runoff evolution and driving factors of Weihe River in Gansu in recent 60 years. Research of Soil and Water Conservation, 30: 183-189. https://doi.org/10.13869/j.cnki.rswc.2023.01.025

[33] Guo, X.F., Ou, T.G., Ma, W.G., Wu, L.B., Zhao, Y.F., Liu, J., Xu,C.Y. (2022). A new random noise attenuation method of seismic signal based on CEEMDAN and wavelet transform. Journal of Geodesy and Geodynamics, 42: 1202-1206. https://doi.org/10.14075/j.jgg.2022.11.019

[34] Hochreiter, S., Schmidhuber, J. (1997). Long short-term memory. Neural Computation, 9(8): 1735-1780. https://doi.org/10.1162/neco.1997.9.8.1735