Bidirectional Segmented Detection of Land Use Change Based on Object-Level Multivariate Time Series

Abstract

High-precision information regarding the location, time, and type of land use change is integral to understanding global changes. Time series (TS) analysis of remote sensing images is a powerful method for land use change detection. To address the complexity of sample selection and the salt-and-pepper noise of pixels, we propose a bidirectional segmented detection (BSD) method based on object-level, multivariate TS, that detects the type and time of land use change from Landsat images. In the proposed method, based on the multiresolution segmentation of objects, three dimensions of object-level TS are constructed using the median of the following indices: the normalized difference vegetation index (NDVI), the normalized difference built index (NDBI), and the modified normalized difference water index (MNDWI). Then, BSD with forward and backward detection is performed on the segmented objects to identify the types and times of land use change. Experimental results indicate that the proposed BSD method effectively detects the type and time of land use change with an overall accuracy of 90.49% and a Kappa coefficient of 0.86. It was also observed that the median value of a segmented object is more representative than the commonly used mean value. In addition, compared with traditional methods such as LandTrendr, the proposed method is competitive in terms of time efficiency and accuracy. Thus, the BSD method can promote efficient and accurate land use change detection.

1. Introduction

Land use patterns have changed significantly in recent decades owing to global changes [1,2]. Although various applications have been used to monitor land use change, many challenges exist in extracting dynamic information on land use change over large areas and with a high frequency [3]. With the development of new remote sensing platforms and sensors, significant progress has recently been made in overcoming technical barriers [4]. Rich historical archives of remote sensing images have also made the long-term detection and modeling of land use change possible [5]. Many studies have been conducted to obtain accurate information on the location, time, and type of land use change using remote sensing image time series (TS) [6,7]. TS analysis has unique advantages, especially with regard to its ability to detect changes in vegetation phenology and land use [8,9].

Among the various TS analysis methods, the pixel-level TS analysis of remote sensing images is the most popular. In this method, a TS is constructed for various phases of every pixel based on a particular index, such as the normalized difference vegetation index (NDVI). Then, by comparing the level of similarity between the TS of the various pixels and that of known change types, it is possible to determine the type and time of change from the pixels [10]. This method has been successfully applied to the TS of various images, including Synthetic Aperture Radar (SAR) [11], Landsat [12], and Moderate Resolution Imaging Spectroradiometer (MODIS) [13]. Dynamic time warping (DTW) is a similarity measure that exploits the temporal distortions between two TSs. The concept of DTW was introduced by Sakoe and Chiba [14] to implement automatic speech recognition. DTW, which measures the optimal global alignment between two TSs, is one of the most frequently used methods to quantify the similarity between two series that are not necessarily equal in size [15]. Research has shown that this algorithm, which is based on nonlinear bending, can achieve high matching accuracy [16]. To replace the Euclidean distance and improve classification accuracy, the DTW distance can be integrated into many classification techniques, including nearest neighbors classifiers, support vector machines, and neural networks [17,18]. In addition, DTW can effectively overcome the time alignment problems caused by missing or low-quality data in the long TS classification process [19]. It has been successfully used in satellite image TS analysis [20], classification [21,22,23], and land-use/land-cover (LULC) mapping [24,25]. Many of the aforementioned methods use a single index to measure the similarity of TS, which highly depends on the accuracy of the selected sample points, and cannot detect the change time. However, given the variety of land use types and the uncertainty of temporal nodes, it is difficult to accurately select samples for all types and times of change for the long-term study of land use change. This, in turn, affects the accuracy of change detection. Moreover, detecting land use changes based on pixel-level TS has stringent requirements for the registration accuracy and radiation correction of the remote sensing images. The results also contain significant salt-and-pepper noise, which limits the application potential of the approach [26,27].

Object-level TS of remote sensing images—another type of analysis method—treats each homogeneous object generated via image segmentation as a whole; this better addresses the aforementioned limitations of pixel-level TS methods [28,29]. Two typical methods are often used for object-level change detection: change detection based on synchronous segmentation and change detection after classification. In the first method, multitemporal data are successively superimposed to generate homogeneous objects of similar shapes, sizes, and positions from multitemporal images. This occurs simultaneously with the segmentation and extraction of image objects. The accuracy of this method is highly dependent on the registration accuracy of the multitemporal data [30]. Additionally, the segmentation results often contain objects that are oversegmented or have incomplete boundaries. This is due to the characteristic parameters of multitemporal images that are simultaneously used during segmentation, increasing heterogeneity. Therefore, postprocessing is required to improve the results. In the second method, object-level classifications of images from different phases are performed separately before a comparative analysis. Then, information on categories, geometries, and spatial contexts of the objects is used to determine the changed objects [30,31].

Many studies focus on analyzing the characteristics of the TS of each segmented object to detect changes in land use [32]. The object-level TS method has advantages over the pixel-level TS method in terms of precision [33,34,35]. Nevertheless, many challenges remain. First, the mean value of all pixels within each object is generally selected to represent the object’s representative feature in the TS. The use of the mean weakens the characteristics of the majority of pixels because pixels close to the object’s edge are less representative. This, in turn, reduces the discriminability between the object and other objects [36,37]. Second, the types and times of land use change vary in general. To ensure accuracy, samples for all types and times of change must be selected and used in the detection of change. This is a difficult task to implement. Therefore, detection methods based on incomplete samples cannot accurately identify the type and time of complex land use change [38,39].

To address the challenges outlined above and detect the types and times of land use change from Landsat images, we propose bidirectional segmented detection (BSD) based on object-level multivariate TS. The proposed method uses the median of an object because the median value of a segmented object is more representative than the commonly used mean value. Based on the three-dimensional (3D) DTW TS similarity measurement method with forward and backward detection, our approach adapts to different types and times of land use change and reduces the difficulty of sample selection. We mainly solve the problem of detecting the change time and change type of land use at the same time without using the changed samples. The main contributions of this study are (1) use of the median instead of the mean as the representative feature of objects and (2) combining forward and backward detection processes to detect different types and times of land use change.

2. Study Area and Data

2.1. Study Area

The study area was the Xinbei District of Changzhou City (31°48′–32°03′N and 119°46′–120°01′E) located in the middle of the Yangtze River Delta, China (Figure 1). It is among China’s first batch of national new- and high-tech industrial development zones. The gross domestic product (GDP) in 2016 was 115.503 billion Chinese Yuan, of which 52% was from secondary industries. Rapid urbanization and industrialization in recent years have led to significant changes in regional land use. Land use changes in the Xinbei District are both typical and diverse. Conversions occur not only between the various subtypes of agricultural lands (e.g., conversion from paddy fields to garden lands) but also between construction lands and agricultural lands or between construction lands and water areas.

2.2. Data

All remote sensing images used in this study were downloaded from the official website of the United States Geological Survey (https://earthexplorer.usgs.gov/). These images were acquired by the Landsat 7 Enhanced Thematic Mapper Plus (ETM+) and Landsat 8 Operational Land Imager (OLI) sensors in 2015–2017 from orbit path 119 and row 38. The images of these years were selected as experimental data because the research area developed rapidly during this period; in this dataset, land conversion was frequent, and the number of images in the study area is sufficient, with few images affected by clouds. Among the images, 21 Landsat images of the study area were selected as experimental data (

Table 1. Landsat image data specifications.

Image Date (Day/Month/Year)	Image Number	Image Type	Cloud (%)	Image Date (Day/Month/Year)	Image Number	Image Type	Cloud (%)
07/02/2015	1	LE7	20.72	26/02/2016	12	LE7	12.78
12/04/2015	2	LE7	22.61	30/04/2016	13	LE7	0.01
14/05/2015	3	LE7	4.76	16/05/2016	14	LE7	0.22
02/08/2015	4	LE7	0.19	27/07/2016	15	LC8	5.05
13/10/2015	5	LC8	4.37	28/08/2016	16	LC8	4..10
08/12/2015	6	LE7	1.09	13/09/2016	17	LC8	24.70
16/12/2015	7	LC8	6.60	02/12/2016	18	LC8	13.19
01/01/2016	8	LC8	3.23	12/02/2017	19	LE7	12.19
09/01/2016	9	LE7	0.77	28/02/2017	20	LE7	6.56
25/01/2016	10	LE7	5.57	08/03/2017	21	LC8	7.66
18/02/2016	11	LC8	9.71

Note: LE7 = Landsat Enhanced Thematic Mapper Plus (ETM+) sensor; LC8 = Operational Land Imager (OLI) sensor.

).

3. Methodology

3.1. Overall Method Concept

We proposed a method based on the analysis of object-level multivariate TS for remote sensing detection in the context of land use change. This method first constructs the object-level NDVI, the normalized difference built index (NDBI), and the modified normalized difference water index (MNDWI) TS. Then, the 3D DTW distance is used to measure the similarity and detect the time and type of land use change via segmentation (Figure 2).

The main process involves the following steps:

(1)

Preprocess remote sensing images: This includes radiation and atmospheric corrections, gap-filling, and image clipping.

(2)

Construct object-level multivariate TS: Remote sensing images of typical phases are selected for the multiresolution segmentation of the images into several homogeneous objects. The NDVI, NDBI, and MNDWI are then used to construct the TS. The median is used instead of the mean as the representative feature of the objects.

(3)

Perform BSD on the segmented objects: The TS of unchanged objects of various land uses are chosen as the sample TS. Then, a sample TS curve most similar to the TS pending detection is screened. Next, forward and backward detections are used to determine the time and type of land use change.

3.2. Preprocessing of Remote Sensing Images

Preprocessing of the images involves radiometric calibration, atmospheric correction, gap-filling, and image clipping. To ensure consistency among different sensors and dates, we performed atmospheric correction and radiometric calibration using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) for each Landsat image [40,41]. Failure of the scan line corrector (SLC) on the Landsat 7 ETM+ resulted in stripes in images acquired after June 2003, which affects the results of the TS analysis. The stripes were filled using the ENVI gap-fill plug-in [42,43,44,45].

3.3. TS Construction Based on Object Characteristics

3.3.1. Multiresolution Segmentation Using TS Images

Prior to TS construction, objects with greater homogeneity should be extracted using the segmentation algorithm from the series of remote sensing images. The segmentation algorithm is a bottom-up growth algorithm [46]. It begins searching for the growth point of an object from the pixel of an image and follows the principle of minimum heterogeneity to find neighboring objects. Before searching for the next object, areas or pixels with the least amount of heterogeneity are merged into one area. Growth is halted upon reaching the maximum heterogeneity value set by the user. When segmenting the series of remote sensing images, multitemporal images are first stacked in time order to form a multiband image; then, the multiband image is segmented.

For multiresolution segmentation, numerous studies have demonstrated the importance of the scale parameter. The scale parameter controls the dimension and size of segmented objects, which may directly affect subsequent results [47,48,49]. In numerous applied studies, land-cover extraction mainly relied on a trial-and-error approach, with segmentation scale parameters determined based on previous experience [50]. However, this approach has been deemed inadvisable [51]; other researchers have since proposed methods to determine optimal segmentation scale parameters [51,52,53]. As one such successful approach to scale optimization, Local Variance (LV) and Rates of Change of LV (ROC-LV) can be combined to determine appropriate segmentation scales [54]; the corresponding Estimate Scale Parameter (ESP) tool was made public for optimizing scale parameters and has been successfully implemented in the eCognition software.

Here, image segmentation was performed using the eCognition software, with which the ESP was used as an evaluation tool to obtain the optimal segmentation scale. The ESP was used to calculate the LV of homogeneous objects in images under different segmentation scale parameters and then to determine the optimal scale parameter for object segmentation using the change rates of the LV. When the change rate was at its largest, the corresponding scale was considered the optimal segmentation scale

The images of all time phases are stacked in chronological order to form a new image with 126 bands, in which the bands selected by the LET7 image are Bands 1–5 and Band 7, and the bands selected by LC8 are Bands 2–7. In the eCognition software, the multiresolution segmentation method was adopted. First, the covariance and variance matrixes were calculated for all pixels of the image according to the band, and then, the correlation coefficient between the bands was obtained [55]. The image layers were set to one for all bands to avoid any bias [56,57,58,59]. Then, the segmentation parameter was selected and the shape parameters and compactness parameters were determined via trial-and-error [58,60,61]. The specific process included first setting the range of shape and compactness to 0.1–0.9 and the change step to 0.1. Then, some buildings and paddy fields in the image were randomly selected and vectorized as the reference segments. Subsequently, various combinations were made, and their results were compared to those of the reference segments. Finally, the ESP tool was used to test the images of the study area and estimate the optimal segment parameters.

3.3.2. Construction of Multivariate TS

When constructing the TS, the selected indices must be capable of clearly distinguishing between various land uses. The five land use types considered in this study were water area, woodland, paddy field, construction land, and dry land (

Table 2. Description of land-cover types.

Land Use/land Covers	Description
Water area	Rivers, lakes, canals, channels, and reservoirs
Woodland	Forests, orchards, and shrubbery
Paddy field	Land used regularly to store aquatic crops such as rice, with obvious seasonal variation characteristics
Construction land	Cities, rural residential areas, roads, railways, etc.
Dryland	Land relying on natural precipitation to cultivate dry crops, along with idle wasteland, bare land, etc.

). NDVI, NDBI, and MNDWI are commonly used indices for analyzing Landsat images [7,60]. NDVI can be used to separate land with dense vegetation cover from that with other uses. NDBI can be used to distinguish construction land. MNDWI is effective at determining water areas and distinguishing them from shadows. Therefore, these three indices were selected for constructing the multivariate TS in this study. The NDVI, NDBI, and MNDWI were calculated using Equations (1)–(3), respectively.

$$\mathit{NDVI}=\frac{\mathit{NIR}-R}{\mathit{NIR}+R}$$

(1)

$$\mathit{NDBI}=\frac{\mathit{MIR}-\mathit{NIR}}{\mathit{MIR}+\mathit{NIR}}$$

(2)

$$\mathit{MNDWI}=\frac{\mathit{Green}-\mathit{MIR}}{\mathit{Green}+\mathit{MIR}}$$

(3)

where NIR, R, and MIR represent the near-infrared red, red, and mid-infrared bands, respectively; and Green is the green band.

Objects generated via image segmentation often contain many pixels with different index values. Such differences are especially prominent in the index values of pixels at an object’s edge versus that of internal pixels. The segmentation results show that there are many mixed pixels in the segmented objects in the construction land. We choose four mixed samples of construction land to compare the mean and median NDVI value of pixels in the objects. Meanwhile, we choose the samples of nonconstruction area (woodland land, water area, dry land, and paddy field) for comparison (

Table 3. Statistical results of mean and median of objects. The red line represents the mixed object (segmentation result), and the yellow line represents the comparative homogenous figure boundary selected manually.

Land-Use/Land-Cover Type	High-Resolution Images of Sample	NDVI of Sample	Statistical Results of Sample
Construction land
Construction land
Construction land
Construction land
Woodland
Water area
Dryland
Paddy field

). For each sample, we manually selected a homogenous object of majority pixels within the segmented boundary (indicated by a yellow line). The experimental results show that the mean and median NDVI values of nonconstruction land are almost the same because the mixed pixels in the segmented objects of nonconstruction land are fewer. For each sample of construction land, the mean NDVI values of a multiresolution segmented object and a homogenous object are quite different, whereas their median NDVI values are almost the same. To further validate the results, we selected 300 segmented objects and corresponding homogenous objects of construction land for statistics. Figure 3 shows that the median NDVI values of segmented objects are closer to the mean NDVI values of corresponding homogenous objects. Therefore, the median value of the segmented object is more representative than the commonly used mean value.

Most existing studies on object-level TS use the mean index of all pixels in an object as the representative feature [61,62]. However, this does not reflect the major characteristics of the object’s index. The median index of pixels is more representative because it better reflects the index characteristics of the majority of pixels in the object and reduces the influence of outliers or noise, especially for construction land (Table 3). Therefore, the median index was used to construct the TS.

3.4. BSD of Segmented Object

3.4.1. BSD Concept

When an object undergoes a land use change, it is necessary to detect the time and type of that change. For shorter time spans (i.e., within 10 years), there is usually one change in land use, such as the conversion of paddy fields to construction land [63,64,65]. For such cases, the land use at the beginning and the end of the TS clearly indicate the type of change. The BSD method proposed in this study can simultaneously detect information in two dimensions: change time and change type.

Figure 4 shows the concepts behind the BSD method. In this example, the land use of the object pending detection changed from Type A to B. In this process, its NDVI also changed, with the NDVI TS containing 21 phases. For Phases 1–9, the NDVI TS curve of the object was similar to Type A, then trended toward Type B for Phases 9–10, and became similar to Type B after Phase 10 without any reversal to Type A. The changing trend of the curve indicates that the object underwent a transition in land use from Type A to B. During the process of change, the rate of change in similarity, which was measured using the DTW distance, increased sharply. It is assumed that the phase of the most abrupt change corresponds to the point at which change occurred. Then, the TS prior to this change, the pre-change TS segment, was extracted. This segment was compared to the TS curve of each of the unchanged samples during the corresponding phase interval. During the comparison, the 3D DTW distance was used to determine similarity [12]. Thus, the type of land use before the change was determined in terms of the type of sample with the highest similarity. Similarly, the new types of land use were determined by comparing the TS curve after change to that of samples with unchanged land uses.

3.4.2. Change Detection Using the BSD Method

The BSD steps are further outlined as follows:

a. Screening of similar land use types

Assume that there are M curves to be detected, N sample curves, and that the series length for both is S. For each TS curve to be detected, m_i (i = 1, 2, …, M) comparisons are made between every subset of TS, starting from the first phase and ending in the middle or end phase, and the corresponding interval of each sample curve. The DTW distance for each comparison is then separately calculated and recorded. The cumulative distance D_n (n = 1, 2, …, N) is obtained as shown in Equation (4).

$$D_n={\left({\displaystyle \sum _{i=3}^S}{d}_{1,i}\right)}_{\mathit{NDVI}}+{\left({\displaystyle \sum _{i=3}^S}{d}_{1,i}\right)}_{\mathit{NDBI}}+{\left({\displaystyle \sum _{i=3}^S}{d}_{1,i}\right)}_{\mathit{MNDWI}}$$

(4)

where i represents the ith phase, which denotes the end of the TS subset, S is the last phase, and $d_{1,i}$ is the distance between the two subsequences corresponding to the first time phase and the ith time phase.

Then, the kth (1 ≤ k ≤ N) sample curve for which D_k is minimum can be identified to be most similar to the TS curve pending detection. The two land uses involved in the change, which are yet to be determined, can be identified via forward and backward detection, as described below.

b. Forward detection

The following equation is used to calculate ${\mathit{\text{r\_D}}}_k$, which is the rate of change of the DTW distance:

$${\mathit{\text{r\_D}}}_k{=[r}_i,r_{i+1}\dots r_{i+j}\dots r_S]$$

(5)

where $r_{i+j}=(d_{1,i+j+1}-d_{1,i+j})/(d_{1,i+j}\left)\right(1\le j\le S-i)$. The maximum value of ${\mathit{\text{r\_D}}}_k$ corresponds to the phase or point of abrupt change, which is the time at which the change in land use begins. This phase is considered the segmentation point, and the curve of the TS prior to this point is extracted as the pre-change TS segment. This segment is then compared to the corresponding TS segment of the samples to determine the type of land use before change.

c. Backward detection

Backward detection starts from the end of the TS; the search is performed in reverse toward the segmentation point detected during forward detection. Based on the rate of change in the DTW distance, the phase or point in time at which the change in land use occurs is detected. The subset of the TS between the end of the phase of change and the last phase is extracted as the TS segment after change. Then, this segment is compared to the corresponding TS segment of the samples to identify the new type of land use.

d. Accuracy evaluation

The aforementioned BSD method was used to detect the types and times of land use change. Then, the detection results were converted into vector data for precision verification. Validation points were selected to calculate the confusion matrix and evaluate the accuracy of the detection results.

Because of the low time resolution of the land use data and Google Earth images, accurate time verification could not be achieved. Therefore, this study reclassified the change time as four layers: unchanged, 2015, 2016, and 2017. By using the random sampling method [66], samples were extracted from the changed layer and from the unchanged layer. As our evaluation focused on the accuracy evaluation of a binary classification problem (changed or unchanged), this sample selection strategy could effectively reduce or even eliminate the impact of errors caused by ground reference data [67]. The two phases of land use data were intersected and analyzed to generate validation layers. The sample points were verified by combining the validation layer with a high-spatial-resolution Google Earth image.

4. Results

4.1. TS Characteristics of Various Land Use Types

Landsat satellite images of the Xinbei District (2015–2017) and two phases of the area’s existing land use maps (2015, 2016) were used to select the TS curves of areas with unchanged land use. First, the land use maps of the two phases were superimposed for analysis. Next, Google Earth historical images were referred to and comparisons were made with the boundaries of the object after multiresolution segmentation. The sample objects for which land use had not changed were selected. Random sampling was adopted to extract 50 samples objects from each of the five types of land use (water area, woodland, paddy field, construction land, and dry land) and the mean value was calculated to form the final sample sequence data. Lastly, the NDVI, NDBI, and MNDWI TS of the sample objects of the land use types were constructed (Figure 5).

Of these indices, MNDWI was the best for distinguishing the different types of land use, followed by NDVI, and then NDBI. For all three TSs, the water areas were highly distinguishable from the other types of land use. Paddy fields and woodlands appeared to be similar in the NDVI and NDBI TSs, but they could be better distinguished in the MNDWI TS. The distinction between paddy fields and construction land in the MNDWI TS was not obvious, but clearer patterns of periodic change for paddy fields were evident in the NDVI TS. Therefore, the NDVI TS could be used to effectively distinguish paddy fields from construction lands. Moreover, the partial segments of the TS for paddy fields, woodland, and construction lands were highly similar. Thus, sufficient phases and the high time resolution of the TS played important roles in accurately distinguishing between various land uses.

4.2. Accuracy Evaluation of Land Use Change Detection

The segment parameters were as follows: maximum-scale = 100, shape = 0.2, compactness = 0.5, and weight of multispectral band = 1. Figure 6 shows the resultant ESP curve.

As small-scale (scale < 30) segmentation can cause object fragmentation in Landsat images, it was not considered in this study. When the segmentation scales were 44, 82, and 95, the local variance of the image’s object reached the local maxima. The curve of the change rate showed that the local variance began decreasing after the local maxima. This result indicates that the heterogeneity of the image’s object reached a maximum before decreasing. The segmentation results shown in Figure 7a–d imply that when the segmentation scale is 95, although the edges of many feature objects are intact, the phenomenon of segmented objects is significant; in other words, different ground classes appear in the same object. The segmentation effect with a scale of 82 is better than that with a scale of 95, and segmented objects still exist. In the segmentation effect results with a scale of 44, although some segmented objects exist, most of them are attributed to oversegmentation. As oversegmentation does not have an adverse impact on the detection results, 44 is the optimal scale in terms of both calculation efficiency and object accuracy. Therefore, the optimal segmentation parameters were selected as follows: segmentation scale of 44, shape factor of 0.2, compactness factor of 0.5, and all bands with a weight of one (

Table 4. Segmentation parameters.

Seg. Level	Scale	Wshape	Wcompactness	Wimage Layer	Number of Objects Produced
Level 1	44	0.2	0.5	1	10,639
Level 2	82	0.2	0.5	1	4572
Level 3	95	0.2	0.5	1	3334

).

The time span of detection ranges from 7 February 2015 to 7 March 2017. The time of detection was changed on 2 August 2015, 25 October 2015, and 16 May 2016. On 28 August 2017 and 12 February 2017, 210 samples (70 samples per year) were extracted from the changed layer and another 200 samples, from the unchanged layer (Figure 8).

The detection results calculated using the confusion matrix had an overall accuracy of 90.49% and a Kappa coefficient of 0.86 (

Table 5. Confusion matrix of change time detection.

Detection Results	Reference Data				User Accuracy (%)
	Unchanged	2015	2016	2017
Unchanged	188	5	2	5	94.00
2015	1	62	6	1	88.57
2016	3	5	60	2	85.71
2017	2	3	4	61	87.14
Producer accuracy (%)	96.91	82.67	83.33	88.41

Overall accuracy: 90.49%. Kappa coefficient: 0.86. Total number of reference pixels: 410.

). The user and producer accuracy both exceeded 80%, where the lowest accuracy for producers occurred in 2016 and for users, in 2015. Misdetection primarily occurred in adjacent years, which may be caused by policy planning and engineering implementation not being synchronized. If only the spatial accuracy is considered and after verifying whether the detection results have changed or not considering the change time, the overall accuracy is 94.88% and Kappa coefficient is 0.89. Overall, the accuracy verification confirms that our method can achieve superior results in land use change detection compared to traditional methods.

4.3. Spatiotemporal Characteristics of Land Use Change

The proposed BSD method identified not only the type of change but also the time when change began. The spatial distributions of the types and times of land use change are shown in Figure 9a,b, respectively.

In terms of spatial distribution, there was a certain correlation in the regions of change. There were clusters of land use change at the centers of cities and towns or close to them. The proximity effect was observed in local areas and mainly concentrated in the towns of Chunjiang, Xuejia, and Luoxi. Chunjiang is the largest township by area in the Xinbei District. It is located on the south side of the Yangtze River and has well-developed water and land transportation systems as well as an economic development zone where several large enterprise groups and many promising industrial companies are located. Consequently, there has been rapid land use change in Chunjiang Town, wherein many paddy fields have been converted into construction lands.

Luoxi town is located adjacent to the international airport of Changzhou (Benniu Airport). It has many highways and provincial roads passing through it. It also has the largest inland river harbor in Jiangsu province, Benniu Harbor. The town has convenient air, water, land, and rail transportation facilities. Recently, many of the surrounding paddy fields and drylands have been converted for the expansion of the Benniu Airport. Xuejia town is adjacent to the northern central business district (CBD) of Changzhou, which has led to rapid conversion of paddy fields and drylands into construction land. In general, land use change is more prominent in Xinbei District because of the rapid development of the economy and transportation systems.

The changes were mainly concentrated in the summer and fall (02/08/2015, 13/10/2015, 16/05/2016, 28/08/2016), accounting for 89.66% of the total land use change. There were fewer areas with changes in winter (12/02/2017), which can be attributed to bad weather and frequent holidays (i.e., New Year’s Day and the Spring Festival), which severely affect construction efficiency. Another possible reason for this is that state approval for construction land is usually granted in the second half of the year, and construction often begins after the Spring Festival, which is usually in February. Therefore, most of the change in land use from nonconstruction to construction land occurred in the spring, summer, and fall.

The major type of land use change from February 2015 to March 2017 involved paddy fields being converted to construction land; this accounted for 73.15% of the total area changed. This was because the original land use in the Xinbei District was predominantly paddy fields. In 2016, paddy fields covered 28.62% of the total area. The paddy fields were mostly distributed in western Chunjiang Town, Xinqiao Town, and Menghe Town. However, the paddy fields converted to construction lands were mainly distributed in Luoxi Town, Xuejia Town, and the central and eastern parts of Chunjiang Town. This indicates that economic development was the main cause of land use change. Other than paddy fields, some dry lands and water areas in Xuejia and Chunjiang Town were also converted to construction land. As the patches of dry land in the Xinbei District were small and scattered near the fringes of the cities, they were less affected by the expansion of construction land.

4.4. Comparative Experiments

4.4.1. Comparisons with Results of LandTrendr Algorithm

According to the analysis of detection results, the primary land use change type is characterized by the conversion of other land use types to the construction land use type. To further study the expansion of construction land in the study area, the LandTrendr algorithm is used. The LandTrendr algorithm has provided accurate results when it has been used to detect forest disturbances, farming abandonment in agricultural land, and the disturbance-recovery of open-pit mines [41]. Wang et al. used the LandTrendr algorithm to study urban impervious surface expansion and achieved ideal results [68]. Therefore, the results of the LandTrendr algorithm were compared with those of our method to verify its suitability.

The LandTrendr algorithm can be divided into the following steps. From the TS trajectory parameters calculated by LandTrendr, the ratio of maximum disturbance value, maximum recovery value, disturbance value, and recovery value was selected as the feature. The sample of construction land expansion was obtained by using the land use change data from 2015, 2016, and 2017; the object with a trajectory of construction land expansion was extracted via a support vector machine. Figure 10 shows the extraction results. The same strategy as detailed in Section 4.2 was used to evaluate the accuracy;

Table 6. Confusion matrix of LandTrendr algorithm.

Detection Results	Reference Data				User Accuracy (%)
	Unchanged	2015	2016	2017
Unchanged	172	11	8	9	86.00
2015	2	59	6	3	84.29
2016	3	6	56	5	80.00
2017	2	4	5	59	84.29
Producer accuracy (%)	96.09	73.75	74.67	77.63

Overall accuracy: 84.39%. Kappa coefficient: 0.77. Total number of reference pixels: 410.

shows the resulting confusion matrix.

As shown in Figure 10a, the BSD results are highly consistent with those of LandTrendr, such that the regions with significant changes are accurately extracted. BSD is less effective than LandTrendr in road detection. In terms of detection accuracy, the overall accuracy of LandTrendr is 84.39%; however, the perception of changing years is not as sensitive as in the BSD method, that is, there are more misclassifications in adjacent years.

4.4.2. Comparison with Change Detection Methods Based on TS Similarity

Traditional methods of identifying land use change based on TS similarity are based on the idea of supervised classification. Samples of land use change are selected before classification. Then, the similarity between the TS of objects pending detection and those of the samples is measured individually. With an appropriate distance threshold, if the similarity is less than the threshold, it is assumed that the same type of land use change has occurred for the pending objects and samples. However, it is difficult to accurately determine the time of change because the samples may not cover all the types and times of land use change, which leads to many false detections.

To compare the BSD method with the TS similarity method, a predominant type of land use change—conversion from paddy field to construction land—was detected using both methods. For the TS similarity method, two sample sets were used. The first set was a sample where land use changed on 13 October 2015, with most conversions involving paddy field land use to construction land use. The second set included five samples wherein land use changed on 02 August 2015, 13 October 2015, 16 May 2016, 28 August 2016, and 12 February 2017. The second set covered all change times detected, as outlined in Section 4.4. The similarity threshold used in the experiment was 2.5 (Figure 11).

Many objects where conversion from paddy field to construction land occurred were not detected using the TS similarity method with the first sample set. The second sample set was better than the first; however, it also resulted in many omission errors. The BSD method outperformed the TS similarity method for both the first and the second sample sets. Additionally, unlike in the TS similarity method, which only detected conversions from paddy fields to construction land, the BSD method also detected conversions from water areas and dry land to construction land. The BSD method avoids the issue of sample inadequacy and can thus be used to detect land use change more comprehensively.

4.4.3. Comparison with Change Detection Based on Mean TS Similarity

Most existing studies on object-level TS use the mean index of all pixels in an object as the representative feature. However, the median index of the pixels is more representative because it better reflects the index characteristics of the majority of the pixels in the object and reduces the influence of outliers and noise. Therefore, the median can increase the difference between land types, as shown in Figure 3. Here, we also compared the difference between these two TS approaches (Figure 12a).

Although the overall accuracy of the mean TS reached 83.17%, the producer accuracy is close to 80%, whereas the user accuracy of the changing results is less than 85% (

Table 7. Confusion matrix of mean results.

Detection Results	Reference Data				User Accuracy (%)
	Unchanged	2015	2016	2017
Unchanged	167	10	10	13	83.50
2015	2	61	3	4	87.14
2016	2	6	56	6	80.00
2017	1	5	7	57	81.43
Producer accuracy (%)	97.09	74.39	73.68	71.25

Overall accuracy: 83.17%. Kappa coefficient: 0.76. Total number of reference pixels: 410.

). There are many large plots that were not successfully detected (Figure 12b). We chose a typical sample object. From Figure 12c, the part of the object that actually changed was a building contained in the dashed line; areas I, II, and III included unchanged surface features. This leads to the difference between the mean TS and the median TS curves shown in Figure 12c. It is clear that the median can better characterize change during change detection.

5. Discussion

5.1. Applicability of BSD Method and Discussion of Multiple Change Detection Method

In land use change detection, the commonly used methods cannot simultaneously determine the change time and type. The BSD method proposed herein can effectively address this problem. The TS curves of objects whose land use does not change were used to identify the type and time of change in the segmented objects. This method can be applied to detect land use change for objects with similar types of change but different change times. In this study, both objects A and B underwent the same type of land use change, from paddy fields to construction land (Figure 13). However, the change time was different for each: object A’s change began on 2 August 2015, whereas object B’s change began on 28 August 2016. The proposed BSD method simply and effectively detects both the type and the time of land use change. It does not need samples of land use change for every type and time of change, and only unchanged samples are needed, which simplifies sample selection.

It must be noted that the proposed BSD method cannot be directly applied to areas with a long-term TS or frequent land use changes. However, it can be used if the long-term TS is divided into shorter periods. Figure 14 shows a schematic diagram of the multiple change detection method. First, we need to use the S-G algorithm to filter the long TS, and then five (or more) parts, S1, S2, S3, a, and b, are obtained using the BSD method. Here, a and b are in the change stage and S1 and S2 are respectively the start state and the final state. Next, the change process subsequence (a, b) is removed, and the BSD method is used to detect the change process of the S3 subsequence; therefore, reverse iteration occurs until all subsequences are detected. Finally, land use changes during the various periods are combined to determine the land use change in a long-term TS.

5.2. Experiment with Long TS

To test the applicability of the algorithm for a long TS, this study selects Landsat images from March to May in 2006–2016 (

Table 8. Long TS Landsat image data specifications.

Image Date (Day/Month/Year)	Image Number	Image Type	Cloud (%)	Image Date (Day/Month/Year)	Image Number	Image Type	Cloud
02/03/2006	1	LE7	0.01	03/04/2012	38	LE7	0.09
18/03/2006	2	LE7	5.71	05/05/2012	39	LE7	0.12
21/05/2006	3	LE7	0.00	21/05/2012	40	LE7	31.79
18/09/2006	4	LT5	0.03	29/11/2012	41	LE7	0.11
23/12/2006	5	LT5	16.04	14/04/2013	42	LC8	4.46
08/01/2007	6	LT5	0.01	12/08/2013	43	LE7	0.00
01/02/2007	7	LE7	0.03	16/11/2013	44	LE7	0.09
21/03/2007	8	LE7	0.08	02/12/2013	45	LE7	0.02
29/03/2007	9	LT5	0.03	10/12/2013	46	LC8	0.95
08/05/2007	10	LE7	0.09	03/01/2014	47	LE7	23.23
03/01/2008	11	LE7	0.15	16/03/2014	48	LC8	0.47
12/02/2008	12	LT5	0.02	27/05/2014	49	LE7	22.16
28/02/2008	13	LT5	0.01	30/07/2014	50	LE7	0.10
23/03/2008	14	LE7	0.03	18/10/2014	51	LE7	20.72
24/04/2008	15	LE7	0.05	26/10/2014	52	LC8	0.09
02/05/2008	16	LT5	0.15	03/11/2014	53	LE7	0.01
05/07/2008	17	LT5	0.03	05/12/2014	54	LE7	0.02
13/01/2009	18	LT5	0.02	13/12/2014	55	LC8	14.19
26/03/2009	19	LE7	14.10	29/12/2014	56	LC8	0.61
11/04/2009	20	LE7	3.94	07/02/2015	57	LE7	20.72
25/08/2009	21	LT5	2.45	12/04/2015	58	LE7	22.61
05/11/2009	22	LE7	14.18	14/05/2015	59	LE7	4.76
21/11/2009	23	LE7	0.02	02/08/2015	60	LE7	0.19
23/12/2009	24	LE7	0.00	13/10/2015	61	LC8	4.37
08/01/2010	25	LE7	0.20	08/12/2015	62	LE7	1.09
30/04/2010	26	LE7	0.00	16/12/2015	63	LC8	6.60
24/05/2010	27	LT5	0.00	01/01/2016	64	LC8	3.23
21/09/2010	28	LE7	0.00	09/01/2016	65	LE7	0.77
10/12/2010	29	LE7	0.01	25/01/2016	66	LE7	5.57
18/12/2010	30	LT5	0.02	18/02/2016	67	LC8	9.71
04/02/2011	31	LT5	0.09	26/02/2016	68	LE7	12.78
01/04/2011	32	LE7	3.61	30/04/2016	69	LE7	0.01
09/04/2011	33	LT5	0.00	16/05/2016	70	LE7	0.22
17/04/2011	34	LE7	0.01	27/07/2016	71	LC8	5.05
25/04/2011	35	LT5	0.48	28/08/2016	72	LC8	4..10
19/05/2011	36	LE7	0.22	13/09/2016	73	LC8	24.7
24/09/2011	37	LE7	4.48	02/12/2016	74	LC8	13.19

) and detects the expansion results of construction land in the study area in 2006–2016 using the BSD method after the S-G filter (Figure 15a). By using the confusion matrix to evaluate the accuracy and the stratified sampling method, 500 samples that change every year and 120 that remain unchanged are taken from the results, for a total of 620 sample points (Figure 15b).

Table 9. Confusion matrix.

Detection Results	Reference Data											User Accuracy (%)
	Unchanged	06	07	08	09	10	11	12	13	14	15
changed	115	0	1	1	1	0	1	0	1	0	0	95.83
06	1	43	3	2	0	0	1	0	0	0	0	86.00
07	2	1	44	1	0	1	0	0	1	0	0	88.00
08	1	0	1	46	1	1	0	0	0	0	0	92.00
09	1	0	0	1	47	1	0	0	0	0	0	94.00
10	2	0	1	0	1	45	1	0	0	0	0	90.00
11	1	0	0	0	2	1	45	1	0	0	0	90.00
12	0	0	0	1	0	2	1	44	1	1	0	88.00
13	1	0	0	0	0	1	2	3	42	0	1	84.00
14	1	0	0	0	0	0	2	2	2	42	1	84.00
15	2	0	0	0	0	1	0	0	2	5	40	80.00
Producer accuracy (%)	91.27	95.56	89.80	86.79	90.57	86.54	84.91	88.00	85.71	87.50	95.24

Overall accuracy: 89.19%. Kappa coefficient: 0.88. Total number of reference pixels: 620.

shows the confusion matrix. The overall accuracy is 89.19% and Kappa coefficient is 0.88; this demonstrates that the BSD method can well detect the expansion and change time points of construction land. Among them, five of the 120 unchanged sample points have construction land expansion in different years, whereas only 12 of the sample points have no construction land expansion. In addition, most of the error samples detected at the change time point are one year, but a few are two years. The user accuracy of the change time points in 2006, 2013, and 2015 is low, which may be due to the fact that they are at both ends of the TS, and the subsequences after segmentation are short and cannot be identified accurately.

Therefore, the BSD method proposed in this study can be applied to a long TS, and it can provide new ideas and methods for land use change detection in a large area over a long time. However, for a smaller time range and more precise land change identification, it is still unable to achieve accurate monitoring. In the future, the data and methods will be improved continuously. For example, future work will focus on accurately detecting the changing time points of long-term and multiple changing geographical elements through additional effective algorithms.

5.3. Influence of Gap-Filled Area on Change Detection

The reliability of the pixel value of the gap-filled area directly affects the accuracy of change detection. As mentioned in Section 3.2, many scholars have proved that the gap-filled image has excellent usability. Figure 16a,b show the comparison of the effect of an ETM + image in the study area before and after gap-filling on 30 April 2015. It can be seen that the visual effect after restoration is excellent.

From the TS perspective, this study discussed the impact of the gap-filled ETM + image on land-cover change detection. Figure 16a,b show the images before and after filling, respectively. According to the results of land use change detection, three main sample pixels of land cover were selected (Figure 17a) and their NDVI TS data were counted. As shown in Figure 17b, the gap-filled pixel values of the water land use type are most similar to the normal pixel values for the same type. The effect of the gap-filled paddy field area is not optimal; however, the overall trend is consistent with the normal paddy field pixels. For TS-based research, as long as the trends are consistent, the gap-filled ETM + image does not have a significant impact on the results. The gap-filled value of construction land is slightly higher than the normal value of construction land; however, its sequence trend is highly consistent with that of normal pixels. In conclusion, it is feasible to apply the gap-filled Landsat ETM + images to TS studies.

6. Conclusions

This study proposed bidirectional segmented change detection based on object-level multivariate TS to detect the type and time of land use change from Landsat images. The study area was Xinbei District of Changzhou City. A total of 21 images from 2015–2017 were selected for the experiment. The main conclusions are as follows:

• The proposed BSD method can effectively detect the type and time of an object’s change. Multitemporal images were segmented and the 3D DTW method was adopted to detect the generated object and samples with unchanged land use via the BSD method. The results show that this method can be used to identify objects that have undergone land use changes with an overall accuracy of 90.49% and a Kappa coefficient of 0.86. Moreover, the method has good scalability for a long TS.
• The median index of a segmented object is more representative than the commonly used mean. Although the overall accuracy of the mean TS reached 83.17%, both the omission and the commission were approximately 15%. Objects generated via image segmentation contain clusters of similar pixels. However, differences in the index values between an object’s edge and internal pixels were especially prominent. The median index of the pixels of objects can better extract the main features of object timing and enhance the stability of timing analysis.
• In the Xinbei District from 2015–2017, the main types of land use change were the conversion of paddy fields, dry land, and water areas to construction land; the change from paddy fields to construction land accounted for 73.15% of the total area changed. In terms of spatial distribution, these changes mostly occurred in areas with rapid economic development. These changes were also near transportation hubs and in urban centers, where the demand for expansion was higher. The times of change were mainly concentrated in summer and fall.

Further study is needed on long-term TS with frequent changes in land use to improve the applicability of object-level TS analysis to complex land use change detection.