Land Cover Change Detection in Ulaanbaatar Using the Breaks for Additive Seasonal and Trend Method

Ulaanbaatar, the capital of Mongolia, has expanded rapidly over the past decade. Insufficient authority is in place to address this expansion, and many residential plots have been developed in the peripheral regions of the city. The aim of this study is to estimate changes in land cover within the central part of Ulaanbaatar, which has been affected by anthropogenic disturbances. The breaks for additive seasonal and trend (BFAST) method is a powerful tool for implementing this study because it is able to robustly and automatically derive the timing and locations of land cover changes from spatio-temporal datasets. We applied the BFAST method for the first time to urban expansion analysis, with NDVI time series calculated from MODIS (MOD09A1 product) during the period 2000–2010. The results show that land cover has changed across 22.51% of the study area, and that the change occurs at a later time with increasing distance from the city center. Bi-temporal high-resolution satellite images of a sample area in 2000 and 2008 confirmed that the detection of land cover changes by BFAST corresponds to areas in which residential development is dominant. This study demonstrates that BFAST is an effective method for monitoring urban expansion. In addition, it increases the applicability of NDVI time series. OPEN ACCESS


Introduction
Rapid urbanization is a paramount issue in global development, and must be continually monitored and controlled.The United Nations forecast a rise in the global population from 7 billion in 2011 to more than 8 billion in 2030, with 60 percent of this population concentrated in urban areas [1].On average, urban areas are expanding twice as fast as the population is growing [2,3].Urban areas are economic, social, cultural, and political centers that serve as the hub of regional development [4], but substantial urban development comes at a price, with haphazard expansion affecting environmental sustainability in and around urban areas.This phenomenon threatens biodiversity [5,6], drives land cover changes, potentially resulting in environmental degradation such as the loss of farmland [7,8], and also the loss of natural resources [9,10].From the perspective of human habitation, one of the most crucial problems faced by developing countries is the rapid growth of informal settlements in peripheral regions [11,12].Such informal settlements develop without planning controls, and the basic infrastructure required for living is significantly lacking [12].
As with other capital cities in developing countries, Ulaanbaatar, the capital of Mongolia, has rapidly expanded in parallel with an increase in population, causing critical development issues [13].Mongolia experienced a drastic transition from a central planning system to a market-based economy in the 1990s.Following a change in policy that previously restricted internal migration, many Mongolians moved from rural areas to Ulaanbaatar seeking employment, an income and education, and to seek a better standard of living [14,15].By 2010, Mongolia's population had doubled from 0.58 million in 1992, and 41% of the total population was concentrated in Ulaanbaatar, with overspill into grasslands surrounding the urbanized area [16].Many residential plots have developed rapidly in the peripheral regions, and are called "ger-areas".These areas comprise formal and informal residential plots characterized by newly-built detached houses and the construction of "gers", which are traditional Mongolian dwellings designed for a nomadic lifestyle, built easily from wood, and covered with felt.This type of development has a unique feature, in that most residential plots are surrounded by wooden fences, called "khashaas" (see Figure 1) [17].Migrants typically claim open land, build khashaas at the property boundary, and finally build a ger or a detached house on the enclosed land [18].
The main reason for the unrestricted development of these ger-areas is the clash of two different legislation frameworks for land management: between the current master plan and the land reform policy.The current master plan of Ulaanbaatar was developed in 2002, and it seeks to challenge the legitimacy and permanency of these regions, which are viewed as being available only for temporary land use, while the land reform seeks to give the ger-areas formal, permanent, and legal status [15].The master plan has been criticized for having insufficient authority to implement directives, and its apparent inability to address the rapid population growth and expansion of ger-areas [15].In apparent opposition to the master plan, the land reform policy, introduced in 2003 to accelerate the development of a free market economy, has allowed Mongolians to own land and to hold renewable, inheritable, and tradable long-term ownership rights for the first time in Mongolian history [19,20].The land reform policy stipulates that each household is allowed to privatize and possess one plot of land at no cost until 2018.As a result of this stipulation, residential plots have become prevalent, particularly in the city's fringes; and 60% of the population of Ulaanbaatar now lives in these peripheral areas [13,15,21].In managing the rapid urban expansion of cities in developing countries, any development needs to be planned and properly monitored in order to maintain internal equilibrium through the sustainable management of natural resources [22].Appropriate strategies for managing urban expansion must be identified and effectively employed [2].It is therefore clear that an analysis of urban expansion would assist regional planners and decision-makers in understanding growth patterns, thereby allowing plans to be made that include certain essential infrastructures [23].Detecting and characterizing changes in land cover is a natural first step towards identifying and understanding the drivers and mechanisms for such change [23].
Many techniques have been developed to detect and analyze land cover change and apply the methods to different regions.The most popular techniques are post-classification comparison, principal component analysis (PCA), and wavelet decomposition [23][24][25][26][27][28][29][30][31][32][33][34][35].The post-classification comparison technique is used to compare images in which the land use/cover has been classified between different periods, by using high-resolution satellite images such as those from Landsat and SPOT.This method is widely used [27][28][29][30], however some pre-processing of the images is necessary due to differences between observing sensors, serious effects of atmospheric disturbances, missing data due to clouds and shadows, and correction of inaccurately observed time spans.PCA is a simple and effective technique for enhancing information about change.However, it usually neglects seasonal variation by intentionally summarizing time series data (e.g., yearly composite images), resulting in the loss temporal change information [25,[31][32][33].In comparison, wavelet decomposition, normally using vegetation indices like the normalized difference vegetation index (NDVI), which is the most commonly used measure of vegetation, is a useful basis for the development of a land cover change analysis [26,34,35].However, until recently, one of the main problems associated with the application of this technique was that it required the determination of several thresholds [26,34,35], and it often led to misleading outputs, caused by the different spectral and phenological characteristics of differing land cover types [23].
Therefore, a new approach is proposed, using the Breaks For Additive Seasonal and Trend (BFAST) method, which is a type of wave decomposition technique [23] developed to overcome the limitations of the methods described above.It is much more useful in estimating land cover change than other techniques, because it is able to robustly and automatically derive both the time and location of land cover changes from NDVI time series, without needing to select a reference period, set thresholds, or define a change trajectory [23,36].To apply this method, high-level MODIS products are useful.These are provided freely by NASA and have been atmospherically corrected to yield surface reflectance data.They are gridded into a map projection with temporal periods of a high frequency [37].These datasets have been applied by several pieces of urban expansion research in the analysis and implementation of land cover detection [34,37,38].Although the BFAST method has been shown to be useful, it has not previously been used in an urban expansion analysis.We propose that it could be applied widely and should therefore be made more available, to allow it to be a practical tool for urban planners when monitoring urban expansion.A robust and feasible method is needed to estimate land cover change in Ulaanbaatar in particular, because urban expansion is rapidly accelerating in this developing city.The aim of this study, therefore, is to estimate changes in land cover in Ulaanbaatar using the BFAST method, and to explore the capacity of the BFAST method to be a tool for urban expansion analysis.

Study Area
The Municipality of Ulaanbaatar consists of nine administrative districts, which are sub-divided into 124 sub-districts, called "khoroos".Two districts, Baganuur and Bagakhangai, and two khoroos in the Songinokhanirkhan district and the Bayansurkh district are excluded from the study area because of their remote and disparate locations.Therefore, the study area comprises 122 khoroos in seven districts, with a total area of approximately 756.25 km 2 (see Figure 2).The study area is divided into three main areas: the city center area; the ger-areas; and others [21].The center of the city is located around Sukhbaatar Square, in front of the governmental palace, and has been extensively developed with high-rise offices and apartment buildings.
The ger-areas have developed mainly in suburban areas, and over recent years have rapidly expanded.Today, migrants from rural areas continue to arrive in the city to claim their own land and build gers and khashaas.It is evident that there is a severe lack of essential urban infrastructure in these suburban areas, such as piped water systems, sanitation facilities, paved roads, public transportation, and heating systems [13].Issues that must be urgently addressed include social and spatial inequality, water supply and sanitation, waste management, flood risk reduction, and air pollution [13].Ger-areas are also found in the northern part of the study area, where developments mainly comprise summer camps, known as "zuslans".These areas originally contained summer houses and second homes for city residents, but the area is now being exploited by new residents, and a large number of settlement clusters have now been established [39].Although residential plots are found both in the city center and the ger-areas, those in the city center have experienced a gradual conversion to apartment blocks [18].Other areas of note within the study area include: parts of the forest, the riverside, grasslands, and industrial areas such as electronic plants, factories, and the airport.The Tuul River, one of the largest rivers in Mongolia, runs to the south of the study area, and divides Ulaanbaatar north and south.
Since natural land covers are not as affected by anthropogenic land cover changes, compared to urban surfaces [34], the forests and the Tuul River are masked from the analysis in order to focus on land cover changes caused only by anthropogenic disturbances.

BFAST Method
BFAST integrates the iterative decomposition of time series data into trend, seasonal, and remainder components, with methods for detecting changes (Figure 3) [29].The model is used iteratively to fit a piecewise linear trend and a seasonal model, given by the equation: where is the observed data at time t, is the trend component, is the seasonal component, and is the remainder component [23].It is assumed that is piecewise linear with segment-specific slopes and intercepts on the different segments of 1 ( 0).Thus, there are m breakpoints * ,…, * , such that: where i = 1, . . ., m and we define * 0 and * .represents the piecewise phenological cycle on different 1 ( 0) segments divided by the seasonal breakpoints, # ,…, # ( # 0 and # ), shown as: where the coefficients are , and , , K is the number of harmonic terms, and is the frequency.We employed the harmonic model proposed by Verbesselt et al. [36] which sets K = 3 in Equation ( 3).This model has a robust approach that avoids noise, and parameters that can be easily used to characterize phenological change [36].The frequency is set as f = 46 for annual observations of an 8-day time series in this study.The remainder component is the remaining variation in the data beyond that defined by the seasonal and trend components [36].
Breakpoints in trend and seasonal components are detected iteratively [23,36,40] as follows: (1) breakpoints * ,…, * are estimated using the residuals-based moving sum (MOSUM) test, and are assessed by minimizing Bayesian information criterion (BIC) from the seasonally adjusted data , where is first found by the STL method [41]; (2) , , and are estimated using robust regression based on M-estimations; (3) breakpoints # ,…, # are similarly estimated by MOSUM and BIC from the de-trended data ; (4) revised is estimated based on the M-estimation; (5) the estimation of parameters is iteratively performed until the number and position of breakpoints are unchanging.
All parameters above are determined automatically.BFAST requires only the parameterization of the minimal segment size between potentially detected breaks [42]; in this study, this was set to two years.

Data Pre-Processing and Implementation
The MODIS data product MOD09A1 is used in this study, containing calibrated reflectance for seven spectral bands within the 400 nm to 2,500 nm spectral region, and surface reflectance quality control flags at a resolution of 500 m pixels.Each MOD09A1 pixel contains the best possible observation during an 8-day period, selected on the basis of high observation coverage, low viewing angle, the absence of clouds or cloud shadow, and low aerosol loading [43].We acquired MOD09A1 data for the period 2000-2010 and extracted only the best quality pixels using the surface reflectance quality control flags.The NDVI was calculated by (band 2 − band 1)/(band 1 + band 2) × 10,000 in each dataset, making a total of 506 NDVI time-series images.Seven images that were lacking, due to the pre-observation of MODIS between 1 January and 25 February 2000, were set to NA, as implementation of the BFAST method requires a complete set of data covering the observed period (i.e., 46 images × 11 years = 506 images are needed).In order to smooth the data and suppress disturbances, a Savitzky-Golay filter from the TIMESAT processing scheme was applied [44].This filter interpolates any missing data and fits disturbed values to neighboring values within a temporal trajectory [45].The process is fast and simple to use over a large area, and works well for NDVI filtering [46,47].We implemented BFAST through the ndvits package in the R statistical software [41,48].Interpolated time-series images were clipped along the boundary of the study area, and datasets were then extracted in chronological order by the ExtractFile function within the R-ndvits package.In this study we found a high frequency of breakpoints in the trend components during the period 2000-2010 (e.g., Figure 3).Since the trend component represents gradual changes due to interannual climate variability or land degradation [40,49], its breakpoints do not assist in the detection of land cover changes on urban surfaces well.Instead, breakpoints in the seasonal component do indicate changes in seasonal trend patterns, and therefore we focus only on the breakpoints in these seasonal components, which could represent a change in land cover.

Results
The BFAST method was implemented in the study area to detect the number and timings of breakpoints in the seasonal components.
The frequency of detection of seasonal components changes during the period 2000-2010 (see Figure 4).Changes are detected in 22.51% of the study area, and are mostly seen just once during the 11-year period.Since the forest and Tuul River areas are masked, the changes are detected within areas affected by anthropogenic disturbances.These areas are distributed spatially in locations governed by the extent of the ger-areas, the internal development of the ger-areas, and land degradation due to anthropogenic activities, such as the development of earthen roads.The main concentration of detected change is around the edge of the city center, corresponding to locations in and around the ger-areas.A smaller number of breakpoints are found in the zuslan areas in the northern part of the study area.The temporal periods of the detected seasonal component changes are shown in Figure 5, and the timing of the changes estimated by BFAST are summarized each year to facilitate a comparison, similar to that performed by Verbesselt et al. [29].In the areas around the city center, change was mostly detected during 2004 and 2005, and it is clear that the later a change is detected, the further that area will be from the city center (see Figure 5).Most of the change occurred within the suburban areas between the years 2005 and 2007.
Significant difficulties were encountered in evaluating the performance of change detection methods, and result from an inability to adequately characterize outcome accuracies [34].Therefore, in order to verify that our detected results are in accordance with actual land cover changes, the spatial distribution of detected areas in a sample area were compared with bi-temporal satellite images from IKONOS in 2000, and Quickbird in 2008 (see Figure 6).Apartments are illustrated in the southern portion, ger-areas in the central part, and mountainous areas in the northern part of the sample area.Each residential plot is colored to facilitate the identification of the extent of the ger-areas.An expansion of the ger-areas, shown by the increase of residential plots towards the northern mountainous area, was clearly observed in the bi-temporal high-resolution satellite images.In comparison with the BFAST data, few breakpoints of seasonal component changes are found in the apartment areas, while many are found in the ger-areas, mostly corresponding to areas experiencing internal development of ger-areas, or the conversion of open lands to ger-area residential plots.Although some pixels are false findings of land cover changes, it is evident here that the spatial distribution of detected changes corresponds well to the areas of expansion.

Discussion
Until recently, the expansion of urban areas from the city center outwards to peripheral areas was often confirmed by observing an urban surface with satellite images [7,27,[50][51][52].However, such an approach to understanding changes in land cover caused by urban expansion has been insufficient, because most studies have applied change detection techniques to less-frequently observed time-series images such as those from Landsat and SPOT [6,11,27,28,32,53].These satellites could not capture the ger-areas precisely due to the insufficient spatial resolution.Furthermore, the patterns of change on urban surfaces are also not clear, especially in terms of their temporal characteristics.However, we used images from a MODIS satellite, which has the advantage of performing observations on a regular basis (although the images are of a lower resolution than those from Landsat and SPOT), and by using the break points detected by BFAST in the seasonal components, we were able to identify when and where the land cover changed.
Results from a BFAST calculation over sample areas of open land, an area converted from open land to ger-areas, and an area of apartment buildings, indicate the different wavelet characteristics of each land cover type (Figure 7).Although it was difficult to interpret the differences between each type of land cover using the observed data and the trend component, the characteristics were captured well by the seasonal component.The ger-area (after the breakpoint in Figure 7b) and the apartment area (Figure 7c) have a notably strong cyclic interannual variation with one peak in each year, while the hilly open areas (Figure 7a, and the period before the breakpoint in Figure 7b) have a less obvious periodicity.The breakpoint shown in Figure 7b indicates the structural change of the seasonal component and gives a profound insight into the changes of land cover from the NDVI time series.In this way, use of the NDVI time series, particularly by BFAST, is able to describe land cover changes.
Our results show that land cover changes occurred at the edge of the city center region in Ulaanbaatar (Figure 4), and that the changes have tended to occur at a later time with increasing distance from the city center (Figure 5).Due to the land reform legislation, the expansion of urban areas in Ulaanbaatar has been accelerating since 2003.In the first stage of policy implementation, only lands at that time under unadmitted possession were privatized to families who lived on land in the ger-areas [20], and only those who lived there had the right to possess land.After this first wave of privatization, open land around the city center became attractive to new migrants, and the ger-areas developed rapidly.The first step to claiming private land in open areas was the building of khashaas to clarify the property boundary, as described above, and as a result a vast number have been erected in the ger-areas [15].The expansion of the ger-areas has also occurred to the north and north-west of the city, which are predominantly mountainous areas with steep slopes [39].
Previously, the BFAST method has only ever been applied in areas of rich vegetation, for instance, with the aim of distinguishing plantation land from grassland and detecting forest fires [23,36,42,54,55].We have here identified that the BFAST method is also able to monitor land cover changes caused by urban expansion, however we acknowledge that the following lessons may be learned from further study.
Firstly, in this investigation, we focused only on the seasonal component changes of BFAST, instead of the trend component changes.It is known that the detection of the trend component change in BFAST is much more sensitive than that of the seasonal component change [23], as illustrated by the larger number of breakpoints in the trend component (areas of detected change were found to comprise 39.17% of the study area by this method).However, the trend component includes abrupt and gradual changes in the NDVI time series [23], which are regarded as disturbances and noise against the seasonal component changes [36].Consequently, we were not able to capture any meaningful information relating to land cover changes from the trend component.Secondly, in any investigation, the study period should be of a sufficient length to grasp the target phenomenon.We set the study period to be 11 years, from 2000 to 2010, and were required to exclude the first two years (2000 and 2001) and last two years (2009 and 2010) in order to set the segment size as two years.When conducting an analysis of time-series satellite images, this limited period is still considered to be relatively short for the detection of long-term land cover change [29,49].Urban expansion is still progressing in Ulaanbaatar and monitoring should be continued in order to manage it.Continuous long-term observations on a regular basis and the construction of a spatio-temporal database would be useful in further understanding interannual spatial phenomena [52].
Thirdly, it is essential to precisely distinguish land cover types.Forests and the Tuul River were masked to focus on the urbanized areas profoundly affected by anthropogenic disturbances, although many changes were also found in the masked areas covered by rich vegetation.The Tuul River and the forests have different characteristics in the NDVI time series, and are more highly affected by natural disturbances than human ones.Although it is currently difficult to acquire precise land classification maps in rapidly developing urban areas, validation data for change detection need to be collected, to enable future analysis to provide sufficient documentation of change events (i.e., before and after) [34].Ground-based land cover maps also need to be developed to understand detailed land use and land cover changes.Finally, an accurate assessment of the result has yet to be developed.An assessment of the identified areas in which land cover has changed could be carried out using traditional methodologies derived from a confusion or error matrix [34,56].However, as mentioned above, time-series spatial information of the entire land cover is insufficient for the assessment of land cover changes, especially in terms of temporal changes, and although IKONOS and Quickbird images give fully detailed land cover data, the frequency of observations and their spatial ranges are limited.New remote sensing instruments such as RapidEye, which has a spatial resolution of 5 m and was launched in 2008, would provide great opportunities for this large-scale assessment as it has wider spatial extents compared to IKNOS and Quickbird.
Although these approaches can be considered in further studies, we believe that the use of BFAST in this investigation has provided meaningful results that are applicable in the monitoring of urban expansion over large geographic regions.In addition, the application of BFAST, developed alongside a maximum utilization of MODIS time series, which cover the globe, and is available free of charge, could facilitate easy analysis in other regions.

Conclusions
Reliable information about land cover and land cover changes caused by urban expansion is clearly needed to enable informed planning decisions [2], because unplanned growth of ger-areas and the unprecedented pace of urbanization in Ulaanbaatar have resulted in unemployment, traffic congestion, air pollution, and other negative environmental impacts [13].The findings in this study could contribute to an understanding of the characteristics of urban expansion and consequently the countermeasures that must be taken to halt or minimize future haphazard urban expansion.
The BFAST method was applied for the first time to urban expansion analysis, and was able to estimate the time and location of land cover changes at the fringes of Ulaanbaatar, simultaneously and automatically.The results are verified by comparing bi-temporal high-resolution satellite images, which show the actual land cover changes caused by anthropogenic disturbances, such as the development of residential plots.The areas of detected change are concentrated around the city center, indicating the high influence of anthropogenic effects.BFAST is shown to be a robust and applicable method for the monitoring of urban expansion, when carefully applied and focused only on anthropogenic activities using a long observation period centered on the target phenomenon.

Figure 1 .
Figure 1.Residential plots surrounded by khashaas in the ger-area.

Figure 2 .
Figure 2. Overview of the study area.

Figure 3 .
Figure 3. Fitted seasonal, trend, and remainder components for MODIS time series (2000-2010), generated by the breaks for additive seasonal and trend (BFAST) method.The dashed lines represent the dates of detected breakpoints, together with their confidential intervals (red).

Figure 4 .
Figure 4. Spatial distribution of the frequency of seasonal component changes.

Figure 5 .
Figure 5. Spatial distribution of the timing of seasonal component change occurrences (showing only the first detected changes).

Figure 6 .
Figure 6.Spatial distribution of detected changes in the sample area and the expansion of the ger-area observed by bi-temporal images in 2000 and 2008.

Figure 7 .
Figure 7. Results of applying BFAST to sample sites: (a) open land; (b) area converted from open land to ger-area; and (c) apartment.