Dynamic Monitoring of Forest Land in Fuling District Based on Multi-Source Time Series Remote Sensing Images

Abstract

Time series remote sensing images can be used to monitor the dynamic changes of forest lands. Due to consistent cloud cover and fog, a single sensor typically provides limited data for dynamic monitoring. This problem is solved by combining observations from multiple sensors to form a time series (a satellite image time series). In this paper, the pixel-based multi-source remote sensing image fusion (MulTiFuse) method is applied to combine the Landsat time series and Huanjing-1 A/B (HJ-1 A/B) data in the Fuling district of Chongqing, China. The fusion results are further corrected and improved with spatial features. Dynamic monitoring and analysis of the study area are subsequently performed on the improved time series data using the combination of Mann-Kendall trend detection method and Theil Sen Slope analysis. The monitoring results show that a majority of the forest land (60.08%) has experienced strong growth during the 1999–2013 period. Accuracy assessment indicates that the dynamic monitoring using the fused image time series produces results with relatively high accuracies.

1. Introduction

The forest is an important subsystem in the ecosystem that provides a variety of valuable timber, raw materials and other supplies necessary to human’s production and life and ecological balance. China, despite its vast territory, lacks forests and greenery, which makes it ecologically fragile. Since the construction of the Three Gorges Dam, the ecological environment on both sides of the Three Gorges reservoir has been changing continuously owing to various natural and human factors. Therefore, it is of great significance to use long-term remote sensing time series imagery to dynamically monitor the forest land resources on both sides of the reservoir.

Long-term time series analysis based on remote sensing have been proven extremely useful in monitoring forest resources changes [1,2,3,4]. Currently, free satellite images are becoming increasing available to the public. For example, beginning in 2008, all Landsat data can be obtained at no charge. Data from a new satellite, Sentinel-2, are also free of charge now [5]. As a result, the study on time series-based analysis has become increasingly popular in the past decade. Landsat data, with its long-term data archive and good data quality, is the most commonly used data for long-term surface monitoring [4,6]. However, Landsat data alone is not sufficient for appropriate time series analysis due to some issues. For example, the Landsat satellites revisit the same area every 16 days, which makes the time interval of the image time series too long. In addition, in some areas, such as the tropics, Landsat satellite data may be contaminated by clouds and shadows. These issues make the available image sequences sparse and reduce the accuracy of dynamic monitoring [6,7].

Multi-source data fusion is a feasible solution to the aforementioned problems. A number of methods for fusing the multi-source data [8,9,10,11,12,13] have been introduced in recent years and successfully applied to detect forest change caused by natural or human factors [1,11,14]. Because Sentinel-2 data are free of charge at the global scale and have similar wavelengths to Landsat data, the area-to-point regression kriging (ATPRK)-based image fusion approach is implemented on Landsat 8 Operational Land Imager (OLI) and Sentinel-2 data to coordinate their spatial resolutions for continuous global monitoring [12]. The Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) [8], developed in 2006 for predicting daily surface reflectance with Landsat spatial resolution and Moderate Resolution Imaging Spectroradiometer (MODIS) temporal frequency, has been widely adopted and became the foundation for many other methods [9,15]. Although recent results of the STARFM method suggest new opportunities for producing remotely sensed data with both high spatial and temporal resolutions, such spatio-temporal fusion methods usually require at least one pair of Landsat-MODIS images (or images from other similar instruments) acquired on the same day to guide the downscaling process of MODIS (or other coarse-resolution data) on other days [8]. This remains a challenge in some regions, such as cloudy areas. In addition, there might be significant uncertainties in the process of downscaling from the observed spatial resolution of 500 m to the target resolution of 30 m. As cloud cover limits optical satellite time series observations, a pixel-based Multi-sensor Time-series correlation and Fusion approach (MulTiFuse) [11] that exploits the full observation density of optical and Synthetic Aperture Radar (SAR) time series has been presented to improve tropical forest monitoring. It has provided an opportunity to use the medium resolution optical and SAR satellite imagery in a beneficial way for improved forest monitoring in tropical regions. Although MulTiFuse is originally applied to Landsat Normalized Difference Vegetation Index (NDVI) and PALSAR data, it is not restricted to medium resolution optical and SAR data and can be used to fuse a variety of time series from other sensors [11].

With the development of remote sensing technology, a large number of remote sensing time series data have been collected. In the meantime, the free medium resolution remote sensing image data is increasingly accessible. Huanjing-1 A/B (HJ-1 A/B) are two small satellites launched by China in 2008 for environment and disaster monitoring and forecasting [16]. Each of them carries two charge-coupled device (CCD) multi-spectral sensors with the same design principle. The scan width of a single sensor is 360 km and the spatial resolution is 30 m. The images acquired by the HJ-1 A/B multi-spectral sensor consist of four multi-spectral bands (B1: 0.43–0.52 µm; B2: 0.52–0.6 µm; B3: 0.63–0.69 µm; and B4: 0.76–0.9 µm) which are similar to the Thematic Mapper (TM) sensor bands 1 to 4. The orbits of the HJ-1 A and HJ-1 B satellites are exactly the same, and the revisiting period after networking is only two days [16]. HJ-1 A/B data have proven useful in monitoring environmental and ecological changes [17,18]. The data is currently available for free. As the HJ-1 A/B has a spatial resolution of 30 m and a revisit frequency of two days, it will be a useful complement to other medium resolution remote sensing data.

This paper takes Fuling District of Chongqing as the study area. Because the interval of a single source time series is very long due to the consistent cloud and fog in this region, we propose a multi-source time series remote sensing image fusion method based on MulTiFuse. In this approach, Landsat data and HJ-1 A/B data which has a good spatial resolution are combined to form a time series imagery of shorter time intervals. The fused time series data are then used in the subsequent dynamic monitoring process to detect forest land change in the study area.

2. Materials and Methods

2.1. Study Area

The Fuling district (29°21′N~30°01′N, 106°56′E~107°43′E) is located in the middle of Chongqing (Figure 1), China, at the end of the Three Gorges Reservoir Area.

Due to the favorable environment, the plant species in this area are abundant, among which the subtropical plants play a dominant role. According to the Investigation Report on the Type II Forest Resources in Fuling District in 2003, the forest land in this area covers about 113,800 hectares, accounting for 38.7% of the total land area of 294,200 hectares [19]. By the end of 2014, the forest land has raised to 130,000 hectares, which is 44.10% of the study area [20].

2.2. Data and Methods

2.2.1. Landsat and HJ-1 A/B Time Series Remote Sensing Image Data

Landsat series satellites operate on a near-polar sun-synchronous orbit, which allows the satellites to cover most of the Earth’s surface in a way that they capture each area of the world at a constant local sun time from the same angle. The geometric accuracy of the Landsat images is high and their spectral information is stable. In this study, we download 26 scenes of Landsat imagery (path 127, row 39/40) from USGS GLOVIS portal [21]. These images are mainly acquired during summer (mostly around September) when the vegetation is abundant, with less clouds and good overall image quality. Data collection time ranges from 1999 to 2013 (

Table 1. Sensor and acquisition time of the Landsat images used in the study.

Sensor	Acquisition Time
Landsat 7 ETM+	1999-10-01
Landsat 7 ETM+	2000-09-01
Landsat 7 ETM+	2002-09-07
Landsat 5 TM	2003-09-18
Landsat 7 ETM+	2004-09-12
Landsat 7 ETM+	2005-09-15
Landsat 5 TM	2006-08-09
Landsat 7 ETM+	2007-09-21
Landsat 7 ETM+	2008-09-23
Landsat 7 ETM+	2009-08-25
Landsat 7 ETM+	2010-08-28
Landsat 7 ETM+	2011-06-28
Landsat 7 ETM+	2013-06-17

). It is worth mentioning that since the Fuling District belongs to the mid-subtropical humid monsoon climate and has a long summer, the data in October is also used when there is a lack of data due to cloud and fog interference.

The HJ-1 satellite is a Chinese satellite constellation for the monitoring and forecasting of the environment and disasters [17]. We selected three images (30 September 2009, 16 September 2010 and 17 August 2012) captured by the CCD1 sensor of the HJ-1 B satellite in this study.

In order to ensure data consistency in each time series, the FLAASH atmospheric correction model is used to perform the radiation correction and the relative radiation normalization for the two time series remote sensing images, based on the pseudo invariant feature normalization method. Function of mask (Fmask) method [22] is applied thereafter to detect the cloud and cloud shadow of the Landsat time sequence and the strip of Landsat 7. The strips are repaired based on the multi-image local adaptive regression analysis model. As the HJ-1 A/B remote sensing image selected in this study have good image qualities, the method of threshold detection and manual correction is adopted for the detection of cloud and cloud shadow.

2.2.2. Verification Data

The validation data of this study are mainly multi-sensor remote sensing image data and field photos. Because the study area is large, in order to obtain a reasonable coverage of the validation samples, we collect the reference data from Landsat time series, Gaofen-1 (GF-1) satellite, Mapping Satellite-1, ZY-3, ZY-1 02C, and historical WorldView and Quickbird very high resolution (VHR) imagery from Google Earth during the study period, together with some field photos.

When selecting the sample data, the ground area of the reference point should be greater than a 30 m, which is the spatial resolution of the Landsat image. The multi-phase image is first registered with the same projected coordinate system as the Landsat imagery. It is then superimposed to facilitate the experimental verification. All together, we collected 54 sample areas including 1131 pixels.

2.2.3. Multi-Source Time Series Remote Sensing Image Fusion Method Based on Space-Time Features

The optimal weighted least squares fitting approach based on pixel, Multi-sensor Time-series correlation and Fusion (MulTiFuse) [11], is employed to fuse the multi-source remote sensing image. The fused image is further corrected with spatial features. The corrected image time series is used for subsequent dynamic monitoring with the Mann-Kendall detection method and Theil Sen Slope Analysis. Finally, an accuracy assessment is performed to evaluate the monitoring results. The data processing and analysis procedure is shown in Figure 2.

● Remote Sensing Image Fusion Based on Time Features

When two time series remote sensing images from different sensors are combined, it is assumed that the time frames of the discrete time series extracted from the reference image sequence and the image sequence to be fused are inconsistent. As shown in Figure 3 [11], the acquisition time intervals of time series $X$ and $Y$ are different, and there is an overlapping period. Firstly, linear interpolation is applied to interpolate the time points of $X$ on $Y$, and the interpolated time series $X_{int}$ is obtained. Similarly, $Y_{int}$ is also obtained. Then, in cases where linear interpolation may get singular values and the length of the time series is inconsistent, the amplitudes are introduced and normalized to obtain the weight of each acquisition time point. Finally, the fitting relationship between the interpolation time series $X_{int}$ and $Y_{int}$ is calculated with the weighted least square fitting regression method [23], as shown in the Equation below:

$$Y_{int}=a+b{X}_{int}$$

(1)

where $a$ and $b$ are constants.

● Fusion Image Correction Based on Spatial Features

The above method provides a good idea for multi-source data fusion from the perspective of time characteristics. Based on this, we further uniformly select and modify singular pixels in spatial features of similar pixels to enhance the fusion effect of the multi-source remote sensing image data.

Unsupervised classification refers to the clustering of pixels with similar spectral or spatial features in multispectral images. In this study, the unsupervised classification is uses to cluster image data as the basis to modify singular pixels. The common algorithms for clustering analysis include Iterative Self-Organizing Data Analysis (ISODATA), K-Means, chain method and so on. In ISODATA, the center of each cluster is determined by iteratively calculating the mean value of the sample pixels, which is similar to the K-Means algorithm. However, ISODATA is more advanced than K-Means with better classification results as the algorithm is self-organizing, i.e., being able to adjust the number of clusters according to intermediate results. This study applies ISODATA to perform cluster analysis on time series images to be corrected.

Based on the results of unsupervised classification of multispectral images, the original image corresponding to the image to be corrected is classified, and then the Mahalanobis distance [24] is calculated to determine whether the pixel of the interpolation image belongs to a specific category. The Mahalanobis distance algorithm is:

$$D_{math}^2={\left(X_i-{\overline{X}}_r\right)}^T{{\displaystyle \sum }}^{-1}\left(X_i-{\overline{X}}_r\right)$$

(2)

where $X_i$ indicates the value of pixel $i$ which is to be evaluated, ${\overline{X}}_r$ indicates the average pixel value of category $r$, Σ is the covariance matrix and $D_{math}^2$ reflects the value dissimilarity of the pixels $i$ and the remaining pixels of the category which they belong to.

Because the Mahalanobis distance assumes that the data is multivariate normal, the distance values above approximate follow chi-square distribution. Here, a 90% quantile point is selected as the criterion for evaluating whether a pixel value is abnormal or not. If a pixel value is greater than the quantile point value, it is regarded as abnormal. In this case, the pixel value is replaced with the average pixel value of the category which it belongs to.

● Mann-Kendall and Theil Sen Slope Analysis

The Theil Sen Slope combination Mann-Kendall detection method is an important trend analysis approach [25], which has been increasingly applied to the study of vegetation change [26,27,28]. In the Mann-Kendall trend test, the following assumptions are made: (1) the null hypothesis $H_0$ is: the time series data to be detected $\left(X_1,X_2,\dots ,X_n\right)$ are randomly with no significant trend; (2) the alternative hypothesis $H_1$ is: for all $k$, $j\le n$, and $k\ne j$, the distribution of $X_k$ and $X_j$ has a tendency to increase or decrease. The calculation of test statistics $S$ is shown as:

$$S={\displaystyle \sum }_{k=1}^{n=1}{\displaystyle \sum }_{j=k+1}^n\mathrm{Sgn}\left(X_j-X_k\right)$$

(3)

$$\mathrm{Sgn}\left(X_j-X_k\right)=\{\begin{array}{ll}+1& X_j-X_k>0\\ 0& X_j-X_k=0\\ -1& X_j-X_k<0\end{array}$$

(4)

When $n>10$, $S$ approximately follows standard normal distribution with its mean being 0 and its variance expressed as $V_{ar}\left(S\right)=n\left(n-1\right)\left(2n+5\right)/18$. The standard normal distribution variable $Z$ is the test statistic for the trend test and is computed as:

$$Z=\{\begin{array}{cc}\frac{S-1}{\sqrt{Var\left(S\right)}}& s>0\\ 0& s=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}}& s<0\end{array}$$

(5)

In the bilateral trend test, for a given significance level $a$, if $-Z_{1-a/2}\le Z\le Z_{1-a/2}$, the null hypothesis $H_0$ is considered to be reliable, indicating that there is no apparent change trend in the time series. If $\left|Z\right|>Z_{1-a/2}$, $H_1$ is considered more reliable, implying that the sequence has a significant increase or decrease trend. For the test statistic $Z$, $Z>0$ means that the time series has an increasing trend, while $Z<0$ suggests a decreasing trend. In the significance test, when the confidence is 90%, $Z_{0.95}$ is equal to 1.64; when the confidence is 95%, $Z_{0.975}$ is equal to 1.96.

When $n<10$, a bilateral trend test is performed based on the test statistic $S$. With a given significance level $a$, if $\left|S\right|\le S_{a/2}$, we accept the null hypothesis. It shows that the time series has no statisticalobvious change trend. Otherwise, we reject the null hypothesis and accept the alternative hypothesis, concluding that the time series is significantly increasing or decreasing. For the test statistic $S$, if $S>0$, it indicates that the time series has an upward trend. If $S=0$, it indicates that there is no apparent trend. If $S<0$, it means that the time series has a decreasing trend.

The Theil Sen slope [29,30] is obtained by taking the median value of the ratio of the gray value (or exponent value) of different years to the time difference as shown in Equation (6). The median value can effectively avoid the noise disturbance. There is no significance test for this value.

$$t=\mathrm{Median}\left(\frac{x_j-x_i}{t_j-t_i}\right)$$

(6)

where $1<j<i<N$ ($N$ represents the length of the time series). When $t>0$, the sequence of grayscale (or exponential values) shows an upward trend. When $t<0$, the sequence presents a downward trend. Larger $t$ value implies greater degree of change, and vice versa.

When combining the Mann-Kendall detection method and Theil Sen Slope analysis to analyze the trend of vegetation coverage, with a given confidence of 95% (i.e., significance level of $0.05$, if $t>0$ and $\left|Z\right|\le 1.96$, the time series exhibits an upward trend but not obvious, it can be considered that the vegetation coverage shows a less upward trend. When $t>0$ and $\left|Z\right|>1.96$, the time series presents a significant upward trend, and the vegetation coverage is considered to be a significant upward trend. When $t<0$ and $\left|Z\right|\le 1.96$, the time series appears an upward trend but no obvious trend. It is believed that the vegetation coverage has a less downward trend. When $t<0$ and $\left|Z\right|>1.96$, the time series shows a significant downward trend, and the vegetation coverage can be considered to decrease significantly.

3. Results

The original images are preprocessed using the methods described in Section 2.2.1. The Landsat image sequence is selected as the reference image sequence, and the HJ-1 B image sequence is used as the image sequence to be fused. Then, a fused image sequence is built based on the image fusion method described in Section 2.2.3. Finally, based on Mann-Kendall trend detection method and the Theil Sen slope analysis, the forest land in the study area is dynamically monitored and the accuracy of the results is evaluated.

3.1. Image Fusion Results

The Landsat remote sensing image sequence and the HJ-1 B remote sensing image sequence are fused according to the method discussed in Section 2.2.3. Figure 4 shows the original images of Landsat and HJ-1 B image sequences and the fused image in a certain area of the study area.

Table 2. Statistical values before and after fusion in a sample region.

Data	Band	Minimum	Maximum	Mean Value	Standard Deviation	Root Mean Square Difference (RMSE)
HJ-1 B 20090930	1	95	163	108.576	5.055	0.330
	2	39	83	47.459	3.312	0.530
	3	36	107	49.337	5.775	0.536
	4	32	112	62.556	7.901	0.637
Landsat 7 20090825	1	76	155	98.643	7.185	-
	2	54	146	76.332	8.121	-
	3	40	191	70.826	13.830	-
	4	25	126	63.228	14.684	-
Fused image	1	76.002	172.125	98.251	7.367	0.156
	2	55.489	182.686	75.652	9.456	0.339
	3	36.738	235.146	69.876	13.564	0.188
	4	22.840	128.591	62.035	13.562	0.396

lists the statistical characteristics of the images in the orange squared region before and after the fusion. Combined with visual effects and statistical results, it can be concluded that the fused remote sensing images are more similar to the Landsat.

3.2. Forest Monitoring Results

Figure 5 gives the NDVI values in forest land of Fuling District from 1999 to 2013. As NDVI is the indicator of greenness, we can use its value to predict forest land in the study area. It can be seen that the forest land has kept a generally upward trend during the study period. From 2000 to 2008, the forest increased steadily, and in 2009, it experienced an abrupt growth. The following two events may explain the above phenomena. In 2002, the policy of returning farmland to forests was fully implemented in the Three Gorges reservoir area in order to protect the ecological environment of the reservoir area and prevent soil erosion [31]. In 2008, Chongqing Municipality launched the Chongqing Forest Project [32,33].

Figure 6 displays the geographic distribution of the forest land change in the study area. The forest land range in the study area is based on the interpretation results of the GF-1 satellite image acquired in 2014 and the SPOT data of the Three Gorges Reservoir Area in 2009. According to the statistical results, there are about 121,297 hectares of woodland in growth, about 8969 hectares in decay, and about 451.68 hectares without apparent changes.

Table 3. Percentages of forest land area with different growth trend.

	Strong Recession	Micro-Recession	Unchanged	Weak Growth	Strong Growth
Percentage (%)	2.15	4.72	0.35	31.98	60.08

lists the percentages of forest land area with different of growth trends. It can be seen that the forest with strong growth possesses the highest percentage, accounting for 60.08% of the total forest area.

The remote sensing image trend detection and analysis can not only detect the long-term growth trend of the forest land, but also can set several thresholds and compare them with the calculated statistics $Z$ to obtain the degree of growth of a certain area. Therefore, it is possible to separate the region with a greater degree of change from the less variable area. From the management point of view, policy support can be given priority to the areas with greater reduction of regional growth.

To assess the accuracy of the monitoring results, 54 sample areas are selected from the high spatial resolution remote sensing images, Google Earth, Landsat time series. A total of 1131 pixels are identified to evaluate the accuracy. Since there is no unanimous method to verify the results of the dynamic monitoring based on time series remote sensing images currently, we adopt the overall accuracy, user’s accuracy and producer’s accuracy, which are commonly used indices for accuracy assessment in remote sensing [34,35,36].

Figure 7 and Figure 8 show two sample areas with forest attenuation during the 1999–2013 period. Figure 7 illustrates the woodland decay caused by human production and construction. Figure 8 shows another forest reduction as a result of a forest fire in 2008. From the NDVI trend chart, the forest land is recovered quickly after the fire.

Table 4. Trend detection accuracy.

Class	Positive Change	Unchanged	Negative Change	Total	User’s Accuracy (%)
Positive change	614	32	23	669	91.78
Unchanged	24	116	41	181	64.09
Negative change	17	8	256	281	91.10
Total	655	156	320	1131	-
Producer’s Accuracy (%)	93.74	74.36	80.00	Overall Accuracy (%)	87.18

lists the confusion matrix of the accuracy assessment. The overall accuracy of the dynamic monitoring is 87.18%, implying a generally high accuracy of the monitoring results. Except for the unchanged, both positive change and negative change have high user’s accuracies and producer’s accuracies, reflecting that the forest monitoring performs well in these two categories.

4. Discussion

Due to the interference of weather conditions (such as cloud and fog) and the revisit time of the satellite itself, the time interval of available image sequence in certain regions are considerably long, sometimes even for years. It creates a major challenge to the dynamic monitoring of forests based on long-term remote sensing sequence images in these regions. The results of this paper show that the fusion of Landsat and HJ-1 A/B images has the potential to improve forest monitoring in Fuling, a mountainous area with many cloudy weathers. This may provide ideas for other areas with similar conditions. HJ-1 A/B is currently in extended service. However, the Sentinel-2 launched in 2015, including multi-spectral bands with spatial resolutions of 10 m, 20 m and 60 m, can be used as substitutions for (or supplements to) future time series.

The optical satellite images used in this study are susceptible to cloud cover, which limits the improvement of the time series. When monitoring rapid changes is need, SAR data can be considered as it is not affected by cloud cover.

5. Conclusions

In this paper, Landsat satellite and HJ-1 B satellite sequence images were fused using a temporal and spatial feature-based multi-source time series remote sensing image fusion method. The fused time series was then used to monitor the long-term growth trend of the forest land in Fuling district of Chongqing, China during the 1999–2013 period based on a combination of the Mann-Kendall trend detection and Theil Sen Slope analysis methods. The monitoring results reflect the spatial-temporal evolutionary characteristics of forest growth conditions. A majority of forest land has continued to grow during the study period. A small number of forest areas in the southeast of the study area showed a decreasing trend. Necessary attention and policy should be given to these decreasing areas by the relevant departments. The confusion matrix shows that the dynamic monitoring of the fused time series performs well with high overall, producer’s, and user’s accuracies.

With the availability of more and more remote sensing image data with different resolutions from different sources in the future, remote sensing image analysis based on multi-source time series will become a long-term research topic and is beneficial to many application areas such as environmental monitoring, disaster risk assessment, agricultural prediction and more. Therefore, it is of great significance to study this research topic both for theoretical and practical purposes.